! KP Notes 17 March 2004.
! 1. All "stop <string>" should be replaced with "call ErrorStop(<string>)". This
!  stops and message goes to compile.out when are running "compile > compile.out".
! 2. Function subset. See the question near the start.
! 3. It would be handy to have a subroutine which, for each $prod or
!     $demand function, gets the s: and t: fields  as fields and which says 
!     whether or not each is present.
!    See KP's attempt - routine GetElasticityInfo added into MPSGE module.
! 4. KP added records arg and labelinfo to type symbol. Added these 
!     into idlist only for types sector, commodity, consumer. [It would
!     be easy to do ditto for parameter if that was desirable.]
!     New routine FindLabelFromLine added to assist with this.
! 5. KP added "Messages" structure for error and warning messages.
! 6. Should pick up error (using undeclared commodity in i: or o: field) 
!      in Tom's part of the code?? See the example below.
! $prod:y            
!  o:py                 q:200
!* NB. This error not picked up ** Need to do so in Tom's part of code??
!* Error:  i:p(SECT)            q:DVHOUS(SECT)
!  i:pc(SECT)            q:DVHOUS(SECT)
! 7. Levels value (dollars) associated q: field is ALWAYS market-price*quantity.
!    [This is true in o: lines and in i: lines]
!     The difference between i: lines and o: lines is as follows:
!      In an o: line, the market price includes the taxes. That is,
!       ps = pm - powerofrax (linear equation)
!      In an i: line, the market price does not include the taxes. That is,
!       ps = pm + poweroftax (linear equation)

! 8. [For Ken who is slow on the uptake on this.]
!  Price in i: or o: field is always the market price.
!   Output. 
!	Market price pm is what households and other firms pay (assuming they are paying no taxes) to buy the output.
!       Firm receives the supply or user price ps. The difference is the tax. ps = pm - to (to is power of tax)
!   Input. 
!	Firm buys goods/factors at market price (initially 1).
!       In addition, firm must pay the taxes. This gives the demand price (?cost price) pd. pd = pm + to (to is power of tax)
!	 The demand price pd is what determines the costs for the firm. 
!
! 

! [17/03/04] Next stuff is OLD.
! KP question for Tom
! It seems that for an "s:" field, %ident is not set up - see results from
!  KP's routine GetProdTopVariableElast
! fixed we think.

! KP question for Tom/Laurent
!  What are the rules for figuring out the market clearing equations?

! KP question for Tom/Laurent
!  In Laurent's basic.mps, has p(s) and p(o) and alias (s,o). Can we insist that
!   no alias is used in the MPSGE description?
! When have p(s) and p(o), GEMPACK needs to know that o is subset of s and/or
!   vice versa.
! Must we exclude elements occurring in the .MPS file?
! Answer. Tom will read declarations and all variables (not quantities - parameters)
!  or sets in : fields must have been declared.

! KP question for Tom/Laurent
! How do I put comments in the .MPS file?
! Answer. ! before trailing comment or * in first column

! KP question.
! What are the implied balancing conditions? Should we test them?
! Answer. The coeffs on the 2 sides of each market clearing equation must be equal.
! That is, when set variables equal to 1, get equality.
! Perhaps calc differences and tell user to check they are small.

! Restriction we may need to impose
! $prod may have two outputs - either two o: fields or 
!  one o: ranging over different set from $prod: field.
! Do we ban this?
! We have to figure out the other restrictions we need to impose. 

!!     Routines for constructing and manipulating general equilibrium
!!      models implemented with MPSGE
!!      @version 1.0
!!      @author Thomas Rutherford

module mpsge

        implicit none

!       Records of source data:

!!      Defines a line of MPSGE source code.

        type line
          character (len=2048) :: text		!! Input records are limited to 2048 characters
          integer :: lineno			!! Line number in the input file
	  type (record), pointer :: record	!! Record to which this line belongs
          type (line), pointer :: next		!! Pointer to the next line
        end type line

!!      Defines a list of input lines as they have been read from the input file.

        type linelist
          type (line), pointer :: line
          type (linelist), pointer :: next
        end type linelist

!!      Elements of the symbol table include various MPSGE variables (commodity
!!      prices, sectoral activity level, consumer incomes and auxiliary variables).
!!      The symbol table also includes sets and parameters which are encountered
!!      in the input file.  This list is maintained to provide some context checking
!!      on input.

        type symbol
! Note by KP. a real constant (eg "2.1" as in "q:2.1") is not added to the symbol table.
!  See the use of the logical function "numeric" in Tom's code.
          character (len=31) :: name		!! Symbol identifier
          character (len=9) :: type		!! One of set, parameter, sector, commodity, consumer, auxiliary
          type (linelist), pointer :: lines     !! Points to lines in which this symbol appears

! Note by KP. ** WARNING. arg may be only set for sector, commodity and consumer. Should also set it for parameter **
	  type (domain), pointer :: arg                 !! Points to the arguments (when relevant)  ! added by KP
	  character (len=255) :: labelinfo      !! labelling info (relevant for sector, commodity, consumer) ! added by KP

          type (symbol), pointer :: next        !! Points to next symbol in the global list
          integer :: dim			!! Dimension of the symbol (undefined=-1,scalar=0, vector=1, etc.)
        end type symbol

!!      A domain is simply a list of sets.  This list may include quoted singletons as well as
!!      references to GAMS sets.  Each element of the set is represented by a pointer to the
!!      set in the symbol table

        type domain
          type (symbol), pointer :: s           !! Pointer to this set in the symbol table
          type (domain), pointer :: next        !! Pointer to the next set in this domain (unassociated if no more dimensions)

          integer :: dim			!! Dimension of this domain
        end type domain

!!      An MPSGE input file consists of a series of records.  Each record is composed of
!!      one or more fields.  An input field consist of a label and some field text.

!!      A field in an MPSGE model consists of two primary components: the label and
!!      the field data.  The field data in turn consists of the identifier,
!!      the identifier argument, the spanning sets and the condition clause.

        type field

!! Verbatim text of leading label

          character (len=255)   :: label

!! Verbatim text of following field 

          character (len=255)   :: data

!! Identifier text (data excluding domain and condition)

          character (len=255)   :: identtext

!! Substring text of the identitier argument domain if present

          character (len=255)   :: argtext

!! Substring text of the spanning domain if present (# notation)

          character (len=255)   :: spantext

!! Substring text of the conditional clause if present

          character (len=255)   :: condtn

!! Pointer to field idenfier if the field is readable. 
!! Fields beginning with "(" are skipped.

          type (symbol), pointer :: ident       

!! Pointer to argument domain if it is readable

          type (domain),pointer :: arg

!! Points to span domain if it is present

          type (domain),pointer :: span         

!! Points to source line

          type (line), pointer :: line          

!! Pointer to next field in the data file.
!! The final field in a datafile is unassociated.

          type (field), pointer :: next
                                                
        end type field

!!      This type describes a _record_ in an MPSGE function, possibly spanning multiple lines
!!      if there are continuations (lines "+" in the first column)

        type record
          type (field), pointer  :: firstfield  !! Points to the first field in this record
          type (line), pointer   :: line        !! Points to the input line where this record begins.
          integer                :: nfield      !! Counts the number of fields in this record
          type (record), pointer :: next        
        end type record

!!      A tax coefficient on a production input or output consists of:
!!      - A tax revenue agent (identified in an A: field)
!!      - An ad-valorem tax rate defined by one of the following:
!!        - A T: field representing a fixed and exogenous tax rate
!!        - A N:/M: field pair representing an endogenous tax.
!!      The tax coefficient consists only of pointers to fields in the input file.  Note that 
!!      the A: field may be on a record which appears prior to the tax rate field.
!!      With endogenous taxes, the N: field defines an auxiliary variable which controls the tax
!!      and the M: field defines a multiplier on the tax rate.

        type taxlist

!!      Pointer to A: field which identifes agent who collects this tax revenue:

              type (field), pointer :: a

!!      Pointer to T: field which defines an exogenous ad-valorem tax rate:

              type (field), pointer :: t

!!      Pointer to N: field indicating the auxiliary variable which determines the
!!      ad-valorem tax rate.  For a given tax entry, either T: or N: will be
!!      used, but not both:

              type (field), pointer :: n

!!      Pointer to an M: field defining the multiplier which is to be applied on 
!!      the previously listed endogenous tax.  An M: field must follow immediately
!!      after an N: field.

              type (field), pointer :: m

!!      Next tax applied on the same production input or output.  This is unassociated
!!      if there are no more applicable taxes:

              type (taxlist), pointer :: next
        end type taxlist

!!      This type defines an input or output in a production block.  The domain associated with
!!      a netput consists of the union of sets in the function declaration block and the 
!!      commodity price.  For example, if we have
!!>
!!      $prod:y(i,t)
!!              o:p(j,t)        q:y0(i)
!!              ...
!!<
!!      then the _domain_ of the netput, _d_, is (_i,t,j_), and the _subdomain_ of the netput, _sd_, is (_j_).
!!      

        type netput
                type (field), pointer    :: f           !! Pointer to the I: or O: field
                type (field), pointer    :: qref        !! Pointer to the associated Q: field in the same record
                type (field), pointer    :: pref        !! Pointer to the associated P: field in the same record
                type (domain), pointer   :: d           !! Pointer to the domain of this netput
                type (domain), pointer   :: sd          !! Pointer to the subdomain of the netput
                type (taxlist), pointer  :: t           !! Pointer to the taxes which apply to this netput
		character (len=1024) :: nestcount
		integer :: nestfields
		logical :: output			!! Logical flag for an output coefficient
                integer :: index                        !! Integer index associated with this netput
        end type netput

        type endowment
                type (field), pointer    :: f           !! Pointer to the E: field
                type (field), pointer    :: q           !! Pointer to the Q: field
		type (field), pointer    :: r		!! Pointer to the R: field
                type (domain), pointer   :: d           !! Pointer to the domain of the endowment
                type (domain), pointer   :: sd          !! Pointer to the subdomain of the endowment
		integer :: index
        end type endowment

!!	List of endowment records

	type endowmentlist
		type (endowment), pointer     :: e	
		type (endowmentlist), pointer :: next
	end type endowmentlist


!!      List of netputs references from the associated nest.

        type netputlist
                type (netput), pointer :: np            !! Netput at this location in the list
                type (field),  pointer :: f             !! Nest assignment field
                type (netputlist), pointer :: next      !! Pointer to next element in list (nul if last input)
        end type netputlist

!!      Nest type characterizes a CES nest in a cost function.

        type nest
                character (len=255) :: ident		!! Identifier text for this nest 
                character (len=255) :: parentid		!! Identifier text for the parent nest
                character (len=255) :: value            !! Text field of the elasticity value
                logical :: output                       !! Logical flag identifies type of nest (output or input)
                type (field), pointer :: f              !! Pointer to the input field which defines the nest
                type (nest), pointer :: parent          !! Pointer to parent nest
                type (nestlist), pointer :: nlist       !! Pointer to parent nestlist entry
                type (netputlist), pointer :: elements  !! Pointer to a list of constituent inputs or outputs
                type (nestlist), pointer :: subnests    !! Pointer to a list of constituent nests
                type (nest), pointer :: next            !! Pointer to next nest in current function
                integer :: index                        !! Integer index associated with this nest
		integer :: level			!! Level in the nesting structure (top level = 0)
		logical :: assigned			!! Flag that all the logic for this nest has been written
	end type nest

!!      Linked list of nests used to represent the nested CES substitution structure:

        type nestlist
                type (nest), pointer :: n               !! Pointer to a nest in this list
                type (nestlist), pointer :: next        !! Pointer to next element of the list
        end type nestlist

!! An MPSGE production or demand function

        type function           
          type (record), pointer :: firstrecord !! Pointer to first input record in this function
          integer               :: nrecord      !! Count of number of records in this function
          integer               :: nfield       !! Count of number of fields in this function
          type (nest), pointer  :: inputs       !! Points to the top level input nest    
          type (nest), pointer  :: outputs      !! Points to the top level output nest    
	  type (field), pointer :: demand		!! Points to the demand field
	  type (endowmentlist), pointer :: endowments	!! Points to the endowments declared in this $demand block
          type (function), pointer :: next      !! Pointer in linked list of all functions
        end type function

        type (line), pointer     :: inputfile   !! Pointer to first line in the linked list of input file lines
        type (function), pointer :: functionlist !! Pointer to the first function in linked list
        type (symbol), pointer   :: idlist      !! Pointer to the first symbol
        type (symbol), pointer   :: idlast      !! Pointer to the last symbol

        character (len=32) :: modelname

        integer :: nfun=0               !! Number of functions in the MPSGE file
        integer, parameter :: iomps=10  !! Unit number of the MPSGE file
        integer, parameter :: iogms=11  !! Unit number of the resulting GAMS file
	integer, parameter :: iomge=12	!! Unit number of the resulting MGE file

	logical :: splitlog=.false.	!! Flag for log output of field parsing operation
	logical :: log=.true.		!! Flag for general log output
        logical :: readvars         !! Flag for current input ot declaration records

! [15/03/04] KP Set these true or false depending on what conversion want to do
	logical ConvertToGAMS, ConvertToGEMPACK, ConvertToMGE
! [20/03/04] KP for flexible setting of name of Input file.
	character (len=512) :: InputFileName

! [19/03/04] Structure "Messages" introduced by KP for handling error and warning messages.
! These used in calls to procedures ErrorStop(<string>) and GiveWarning(<string>)
!  and in InternalErrorStop(<string>).
! Note that "\n" inside <string> means to start a new line.
! Idea is that put PreLines before the message <string> and PostLines after it.
! If PreLine(K) is all blank, ignore PreLine(K+1...)
! If PosLine(K) is all blank, ignore PostLine(K+1...)
! Call routine InitializeMessages to set all lines blank
	type MessageStuff
	  character (len=255) :: PreLine(4)
	  character (len=255) :: PostLine(4)
	end type MessageStuff

	type (MessageStuff) :: Messages 

! module used in conversion to GAMS GMS file

	type gmsrecord
	  character (len=1024) :: left, right
	  character (len=32) :: type
	  type (gmsrecord), pointer :: next
	end type gmsrecord

	type (gmsrecord), pointer :: gmsrecords, lastrecord

	character (len=6), parameter :: fmt = '(1x,a)'

contains

!!      Master routine to read an MPSGE input file:
!!        - Open the file
!!        - Read one line at a time into a linked list for listing.
!!        - Locate functions and set up pointers to the constituent records.
!!        - Parse records into fields.
!!
!!	None of the logical structure of the underlying economic model is
!!	formulated at this stage.

        subroutine readmps
        integer :: lineno
	type (line), pointer :: L
        type(function), pointer :: fn
        type(record), pointer :: r
        type(field), pointer :: f, ff
        character (len=2048) :: input=' '
        type(function), pointer :: fn2 ! KP

        nullify(inputfile)
	nullify(functionlist)
	nullify(f)
	nullify(L)

	if(ConvertToGAMS) then   ! KP
	  write(iogms,fmt) '$ontext'
	end if

        lineno = 0
        do

!       Read a line until we reach one which is not a comment:
!	L is nullified at end of file.

	  call readline(L, input, lineno)
	  if (.not.associated(L)) exit

!       Read a new function?

          if (input(1:1).eq.'$') then
            nfun = nfun + 1
	    call InitializeMessages('pre',0)  ! KP
	    Messages%PreLine(1) = 'Processing section "'//trim(input)//'"'  ! KP
	    write(*,*)'NFUN=', nfun, ' INPUT=',trim(input) ! KP

            if (associated(functionlist)) then
	      call newfunction(fn%next); fn => fn%next
	      fn2 => functionlist   ! KP
            else
	      call newfunction(functionlist); fn => functionlist
            endif

!	Set a flag here if we are reading variable declarations to avoid
!	pulling the plug if we can't find a colon:

	    if (eqv(input,'$sectors:') .or. eqv(input,'$commodities:') .or. &
		eqv(input,'$consumers:') .or. eqv(input,'$auxiliary:') .or. &
		eqv(input(1:12), '$constraint:') ) then  ! KP [will be text after "$constraint:", hence (1:12)]
	      readvars = .true.
	    else
	      readvars = .false.
	    end if
          else
            if (.not.associated(functionlist)) &
              call ErrorStop('Error == MPSGE input must begin with a function.')
          end if

!       Read a new record?

          if (input(1:1).ne.'+') then
            if (associated(fn%firstrecord)) then
	      call newrecord(r%next);  r => r%next; 
            else
	      call newrecord(r); fn%firstrecord => r; 
            end if
            r%line => L
            fn%nrecord = fn%nrecord + 1
          end if
	  nullify(f)
	  L%record => r

!	Keep pulling fields off the input record until there are no more:

          do while (input.ne.' ')

	    call newfield(ff)

!	    if (associated(lastfield)) lastfield%next => ff

!	This logic only applies to the first field extracted from the
!	line:

            if (input(1:1).ne.'+') r%firstfield => ff
	    if (associated(f)) ff%next => f

            f => ff
            f%line => L
            r%nfield = r%nfield + 1
	    fn%nfield = fn%nfield + 1

	    call readfield(input,f)
          end do

        end do

	if(ConvertToGAMS) then   ! KP
          write(iogms,fmt) '$offtext'
	end if

	if(log) write(*,*) ' NFUN= ', nfun    ! KP

	call InitializeMessages('pre',0)  ! KP

	return
        end subroutine readmps

!!	Read to the next non-blank line from the input file, incrementing
!!	the line count:

	subroutine readline(L, input,lineno)
	type (line), pointer :: L
        character (len=*) :: input
	integer :: lineno
	integer :: ic

	input = ' '
	do while (input.eq.' ') 
          read(iomps,'(a)',end=100) input
	  lineno = lineno + 1

	  if(ConvertToGAMS) then   ! KP
            write(iogms,'(a,i0,a,t12,a)') '[',lineno,']',trim(input)
	  end if
            
	  if (associated(inputfile)) then
	    call newline(L%next)
	    L => L%next
	  else
	    call newline(inputfile)
	    L => inputfile
	  endif
	  L%text = input
	  L%lineno = lineno
            
	  ic = index(input,'!')
	  if (ic.gt.0) input(ic:) = ' '
	  call tab2blank(input)       ! Convert tabs to blanks:
	  if (input(1:1).eq.'*') input = ' '
	end do
	return
100	nullify(L)
	return
	end subroutine readline

!!	Read the last field from a line of input, using the 
!!	colon to locate the label and the data components:	

	subroutine readfield(input,f)
	type (field), pointer :: f
        character (len=*) :: input
	integer :: ic

	ic = index(input,':',back=.true.)
!	if ((.not.readvars) .and. (ic.eq.0)) stop 'No colon?'
	if ((.not.readvars) .and. (ic.eq.0)) call ErrorStop( 'No colon?')
	f%data = adjustl(input(ic+1:))
	if (ic.eq.0) then
	  f%label = ' '
	else
	  input(ic:) = ' '
	  ic = index(input(1:ic-1),' ',back=.true.)
	  f%label = input(ic+1:)
	  call lower(f%label)
	end if
	input(ic+1:) = ' '
	if (input.eq. '+') input = ' '
	return
	end subroutine readfield

!!      Fold a character string to lower case.

        subroutine lower(word)
        IMPLICIT NONE
        character (len=*) :: word  !! Character to be placed in lower case.
        integer :: i
        do 100 i=1,len(word)
          if (word (i:i).eq.'A') word (i:i) = 'a'
          if (word (i:i).eq.'B') word (i:i) = 'b'
          if (word (i:i).eq.'C') word (i:i) = 'c'
          if (word (i:i).eq.'D') word (i:i) = 'd'
          if (word (i:i).eq.'E') word (i:i) = 'e'
          if (word (i:i).eq.'F') word (i:i) = 'f'
          if (word (i:i).eq.'G') word (i:i) = 'g'
          if (word (i:i).eq.'H') word (i:i) = 'h'
          if (word (i:i).eq.'I') word (i:i) = 'i'
          if (word (i:i).eq.'J') word (i:i) = 'j'
          if (word (i:i).eq.'K') word (i:i) = 'k'
          if (word (i:i).eq.'L') word (i:i) = 'l'
          if (word (i:i).eq.'M') word (i:i) = 'm'
          if (word (i:i).eq.'N') word (i:i) = 'n'
          if (word (i:i).eq.'O') word (i:i) = 'o'
          if (word (i:i).eq.'P') word (i:i) = 'p'
          if (word (i:i).eq.'Q') word (i:i) = 'q'
          if (word (i:i).eq.'R') word (i:i) = 'r'
          if (word (i:i).eq.'S') word (i:i) = 's'
          if (word (i:i).eq.'T') word (i:i) = 't'
          if (word (i:i).eq.'U') word (i:i) = 'u'
          if (word (i:i).eq.'V') word (i:i) = 'v'
          if (word (i:i).eq.'W') word (i:i) = 'w'
          if (word (i:i).eq.'X') word (i:i) = 'x'
          if (word (i:i).eq.'Y') word (i:i) = 'y'
          if (word (i:i).eq.'Z') word (i:i) = 'z'
          if (word (i:i).eq.char(9)) word(i:i) = ' '
100     continue
        return
        end subroutine lower

!!      Fold a character string to upper case.

        subroutine upper(word)
        IMPLICIT NONE
        character (len=*) :: word  !! Character to be placed in upper case.
        integer :: i, iacharpos, newpos
	character :: ch1
        do 100 i=1,len(word)
	  ch1 = word(i:i)
	  iacharpos=IACHAR(ch1)
! "a" is ASCII character 97 and "z" is ASCII character 122
! "A" is ASCII character 95 and "z" is ASCII character 90
	  if(iacharpos.ge.97.and.iacharpos.le.122) then
	    newpos = iacharpos - 32
	    word(i:i) = ACHAR(newpos)
	  end if
          if (word (i:i).eq.char(9)) word(i:i) = ' '
100     continue
        return
        end subroutine upper

!!      Convert tabs on an input line to spaces:

        subroutine tab2blank(line)
        character (len=*) :: line  !! Line to scan for tabs
        character :: tab=char(9)
        integer :: k
        do k=1,len_trim(line)
          if (line(k:k).eq.tab) line(k:k) = ' '
        end do
        return
        end subroutine tab2blank

!!      Evaluate the dimension of a domain by counting
!!      the number of set identifiers in the linked list
!!      defining the domain:

        integer function dimension(d)
        type (domain), pointer :: d     !! Domain to be evaluated
        type (domain), pointer :: dd    
        dimension = 0
        if (.not.associated(d)) return
        dd => d
        do while (associated(dd))
          dimension = dimension + 1
          dd => dd%next
        end do
        end function dimension

!!      Generate an echoprint of the file.  This is currently very preliminary,
!!      and will later be recoded to produce a more carefully formatted report.

        subroutine echoprint

        type(function), pointer :: fn
        type(record), pointer :: r
        type(field), pointer :: f
	integer :: i

        fn => functionlist
        nullify(f);
        do while (associated(fn))
          write(*,fmt) ' ===== '
	  r => fn%firstrecord
	  do while (associated(r))
	    f => r%firstfield
	    do while (associated(f))
	      write(*,fmt) ' '//trim(f%label)//':'//trim(f%data)
	      write(*,fmt) ' Identifier text: '//trim(f%identtext) ! KP
	      write(*,fmt) ' Argument text: '//trim(f%argtext) ! KP
	      f => f%next
	    end do
	    r => r%next
	  end do
	  fn => fn%next
        end do

	return
	end subroutine echoprint

!!      Parse all fields in the file into constituent elements:
!!      - Field identifier text
!!      - Identifier argument text
!!      - Spanning domain text
!!      - Condition text.

        subroutine parse
        type(function), pointer :: fn
        type(record), pointer :: r
        type(field), pointer :: f
        type(symbol), pointer :: s

        fn => functionlist
        nullify(f);
        do while (associated(fn))
	  r => fn%firstrecord
	  do while (associated(r))
	    f => r%firstfield
	    do while (associated(f))

!	Sort out the data:

	      call split(f%data, f%identtext, f%argtext, f%spantext, f%condtn)

! KP (Temporary)
	      write(*,*)'Parse: f%label=', trim(f%label), ' f%identext=', trim(f%identtext)

!       Read domain for each field:

	      if (f%argtext .ne.' ') call readdomain(f%argtext,  f%arg, f%line)
	      if (f%spantext.ne.' ') call readdomain(f%spantext, f%span, f%line)
! Addition 29/01/04 by KP to set f%ident in all cases
	      if (f%identtext.ne.' ') then
		call findid(f%identtext, s)
		f%ident => s

! KP (Temporary?)
!.NO	        if(.not.associated(f%ident)) then
!.NO		  call InitializeMessages('all')
!.NO	          call InternalErrorStop('parse: f%ident is not associated.')
!.NO	        end if
		
	      end if
	      f => f%next
	    end do
	    r => r%next
	  end do
	  fn => fn%next
	end do
	return
        end subroutine parse

!! Return pointer to given identifier

        subroutine findid(name, id)             
        character (len=*), intent(in) :: name   !! Identifier for the symbol to be located
        type (symbol), pointer :: id            !! Pointer into the master list for the located symbol.
        type (symbol), pointer :: s

        nullify(id)
        s => idlist
        do while (associated(s))
          if (s%name==name) then
            id => s
            return
          end if
          s => s%next
        end do
        return 
        end subroutine findid

!! Introduce a new identifier in the master list.

        subroutine addid(name, type, s, L)      

!!      Name of the identifier to be added to the master list

        character (len=*), intent(in) :: name

!!      Description of the type of item to be introduced:
!!      - _set_: a GAMS set identifier
!!      - _parameter_ : a GAMS parameter identifier
!!      - _sector_ : an MPSGE sector (acivity variable)
!!      - _commodity_ : an MPSGE commodity (price variable)
!!      - _consumer_ : an MPSGE consumer (income variable)
!!      - _auxiliary_ : an MPSGE auxiliary variable

        character (len=*), intent(in) :: type
        type (symbol), pointer :: s
        type (line), pointer :: L

        if (.not.associated(idlist)) then
	  call newsymbol(idlist)
          s => idlist
        else
	  call newsymbol(idlast%next)
          s => idlast%next
        end if
        s%name = name
        s%type = type
        allocate(s%lines);
        s%lines%line => L
        nullify(s%lines%next)
        nullify(s%next)

! Additions by KP
	nullify(s%arg)  ! KP
	s%labelinfo= ''  ! KP

        idlast => s
        return
        end subroutine addid

!	Create a copy of a domain.

        subroutine domaincopy(d1,d2)            
        type (domain), pointer:: d1	!! Domain to be copied
        type (domain), pointer :: d2	!! New domain to be created as a copy of d1
        type (domain), pointer :: d,dd
        if (.not.associated(d1)) return

	call newdomain(d2)
        d  => d1; dd => d2
        do while (associated(d))
          dd%s => d%s
          if (associated(d%next)) then
	    call newdomain(dd%next); dd => dd%next
          endif
          d => d%next
        end do
        return
        end subroutine domaincopy

!!	Add an entry in the linked list of line numbers which refer to a given symbol.

        subroutine addlineref(L,s)

        type (symbol), pointer  :: s  !! Symbol which is referenced
        type (line), pointer :: L	!! Line in which the symbol is referenced
        type (linelist), pointer :: LL

        LL => s%lines
        do while (associated(LL))
          if (L%lineno.eq.LL%line%lineno) return
          if (.not. associated(LL%next)) exit
          LL => LL%next
        end do
        allocate(LL%next)
        LL => LL%next
        LL%line => L
        nullify(LL%next)
        return
        end subroutine addlineref

!! Add a set element to the domain

        subroutine domainadd(name,dm, L)     
        character (len=*), intent (in) :: name	!! Identier of set to be added to the domain
        type (domain), pointer :: dm  !! Domain to which this set is to be introduced
        type (line), pointer :: L  !! Line in which the set reference was made.

        type (symbol), pointer  :: s
        type (domain), pointer :: d

!       Find what set we are introducing

        call findid(name, s)
        if (associated(s)) then
          if (s%type.ne.'set') stop 'Error -- not a SET'
          call addlineref(L, s)
        else
          call addid(name,'set',s,L)
        end if

        if (.not. associated(dm)) then
	  call newdomain(dm); dm%s => s; 
          return
        endif
        
!       Add this set at the end of the domain list:

        d => dm
        do 
          if (.not.associated(d%next)) exit
          d => d%next
        end do
	call newdomain(d%next)
	d=>d%next; d%s => s; 

        return
        end subroutine domainadd


!!  Logical function which reports whether a set with a given name
!!  is in a domain:

        logical function indomain(s,dm)
        type (symbol), pointer	:: s	!! Set to be located
        type (domain), pointer	:: dm	!! Domain in which the set is to be located.

        type (domain), pointer :: d

        if (.not. associated(dm)) then
          indomain =.false.
          return
        endif

        d => dm
        do while (associated(d))
          if (associated(d%s,s)) then
            indomain = .true.
            return
          endif
          d => d%next
        end do
        indomain = .false.
        return
        end function indomain

!!   Tets whether one domain (_d1_) is a subst of another (_d2_).

        logical function subset(d1,d2)
! KP really tests if one domain is a subdomain of another domain - that is,
!   sets in d1 are all somewhere amongst the sets in d2, and in the same order
! Is this correct? 
        type (domain), pointer :: d1	!! Subset domain
	type (domain), pointer :: d2	!! Superset domain 
	type (domain), pointer :: d, dd
        d  => d1; 
        do while (associated(d))
          dd => d2;
          subset = .false.
          do while (associated(dd))
            if (d%s%name.eq.dd%s%name) then
              subset = .true.
              exit   ! From "do while" ! KP
            end if
            dd => dd%next
          end do
          if (.not.subset) return
          d => d%next
        end do
        return
        end function subset

        logical function IsSubdomainIgnoreOrder(d1,d2)   ! Added by KP
! KP really tests if one domain is a subdomain of another domain - that is,
!  sets in d1 are all somewhere amongst the sets in d2. Order does not matter.
        type (domain), pointer :: d1	!! Subset domain
	type (domain), pointer :: d2	!! Superset domain 
! Local
	type (domain), pointer :: d, dd
	logical :: DFoundInDD
	
        d  => d1;
	IsSubDomainIgnoreOrder = .false. 
        do while (associated(d))
          dd => d2;
	  DFoundInDD = .false.
          do while (associated(dd))
            if (trim(d%s%name).eq.trim(dd%s%name)) then
	      DFoundInDD = .true.
              exit   ! from "do while"
            end if
            dd => dd%next
          end do
          if (.not.DFoundInDD) then
! temporary
		write(*,*)'IsSubdomainIgnoreOrder false. d%s%name is: ', d%s%name
		IsSubDomainIgnoreOrder = .false. 
		return
	  end if
          d => d%next
        end do
	
! Here result is true
	IsSubDomainIgnoreOrder = .true.
 
        return
        end function IsSubdomainIgnoreOrder

!!      Split up an input field into the constituent elements
!!      - Field identifier text
!!      - Identifier argument text
!!      - Spanning domain text
!!      - Condition text.
!!
!!      For example:
!!>
!!      i:p(i,t)#(s)$beta       
!!<
!!      parses into:
!!>
!!      Field identifier: p
!!      Argument text: i,t
!!      Spanning domain text: s
!!      Condition text: beta
!!<
        subroutine split(field,ident,arg,span,condtn)
        character (len=*), intent(in)  :: field
        character (len=*), intent(out) :: ident,arg,span,condtn
        integer :: ic, np

        ident = ' '
        arg = ' '
        span = ' '
        condtn = ' '

!       Split field text into text fields for the identifier, the argument domain,
!       the span domain and conditional clause:

        np = 0
        ident = field
        do ic=1,len_trim(field)
          if (field(ic:ic).eq.'(') np = np + 1
          if (field(ic:ic).eq.')') np = np - 1
          if (ic.gt.1 .and. field(ic:ic).eq.'$' .and. np.eq.0) then
            condtn = field(ic+1:)
            ident = field(1:ic-1)
            exit
          end if
        end do
        if (np.ne.0) stop ' *** Parethesis mismatch.'

        ic = index(ident,'#')
        if (ic.gt.0) then
          span = adjustl(ident(ic+1:))
          ident(ic:) = ' '
        end if

        ic = index(ident,'(')
        if (ic.gt.1) then
          arg = ident(ic:)
          ident(ic:) = ' '
        end if

	if (splitlog) then
	  write(*,fmt) ' --------------'
	  write(*,fmt) ' field  = '//trim(field)
	  write(*,fmt) ' ident  = '//trim(ident)
	  write(*,fmt) ' arg    = '//trim(arg)
	  write(*,fmt) ' span   = '//trim(span)
	  write(*,fmt) ' condtn = '//trim(condtn)
	end if

        return
        end subroutine split

!!	Read a domain and parse into a series of set elements.  

        subroutine readdomain(dtext,d,L)
        character (len=*), intent(in) :: dtext	!! Domain to be parsed
        type (domain), pointer :: d  !! Pointer to the domain list
        type (line), pointer :: L  !! Pointer to source line
        character (len=255) :: textfield
        character (len=31) :: setid
        integer :: ic

        textfield = dtext

!       Check for leading and closing parentheses (required):

        if (len_trim(textfield).le.2) then
          stop ' *** READDOMAIN: invalid domain text.'
        endif

        if (textfield(1:1).ne.'(') then
          stop ' *** READDOMAIN: opening parenthes in domain text is missing.'
        endif

        if (textfield(len_trim(textfield):len_trim(textfield)).ne.')') then
          stop ' *** READDOMAIN: closing parenthesis in domain text is missing.'
        endif

        textfield = adjustl(textfield(2:len_trim(textfield)-1))

!       Check that there are no extra paretheses:

        if (index(textfield,'(') + index(textfield,')')  + &
            index(textfield,';') + index(textfield,':').gt.0) then
          stop ' *** READDOMAIN: embedded parenthesese or semicolons in domain text.'
        endif

        call lower(textfield)
        do while (textfield.ne.' ')
          ic = index(textfield,',')
          if (ic.eq.0) then
            setid = adjustl(textfield)
            textfield = ' '
          else
            if (ic.eq.1) stop ' *** READDOMAIN: empty set specification.'
            setid = adjustl(textfield(1:ic-1))
            textfield = textfield(ic+1:)
          end if
          call domainadd(setid, d, L)
        end do

        return
        end subroutine readdomain

!!	Allocate a new line and assign default values

	subroutine newline(L)
	type (line), pointer :: L
	nullify(L)
	allocate(L)
	nullify(L%next)
	L%lineno = 0
	L%text = ' '
	return
	end subroutine newline

!!	Allocate a new symbol and assign default values
!!	for all the symbol characteristics

	subroutine newsymbol(s)
	type (symbol), pointer :: s
	nullify(s)
	allocate(s)
	s%name = ' '
	s%type = 'undefined'
	nullify(s%lines)
	nullify(s%next)
	s%dim = -1
	nullify(s%arg)  ! KP
	s%labelinfo = ' '  ! KP
	return
	end subroutine newsymbol

!!	Allocate a new function pointer with default values

	subroutine newfunction(fn)
	type (function), pointer :: fn
	nullify(fn)
	allocate(fn)
	fn%nfield = 0
	fn%nrecord = 0
	nullify(fn%next)
	nullify(fn%inputs)
	nullify(fn%outputs)
        nullify(fn%endowments)
        nullify(fn%demand)
	nullify(fn%firstrecord)
	return
	end subroutine newfunction


!!	Allocate a new domain pointer with default values

	subroutine newdomain(d)
	type (domain), pointer :: d
	nullify(d)
	allocate(d)
	nullify(d%next)
	nullify(d%s)
	d%dim = -1
	return
	end subroutine newdomain

!!	Allocate a new field with default values

	subroutine newfield(f)
	type (field), pointer :: f
	nullify(f)
	allocate(f)
	f%label = ' '
	f%data = ' '
	f%identtext = ' '
	f%argtext = ' '
	f%spantext = ' '
	f%condtn = ' '
	nullify(f%ident)
	nullify(f%arg)
	nullify(f%span)
	nullify(f%line)
	nullify(f%next)
	return
	end subroutine newfield

!!	Allocate a new record with default values

	subroutine newrecord(r)
	type (record), pointer :: r
	nullify(r)
	allocate(r)
	nullify(r%firstfield)
	nullify(r%line)
	r%nfield = 0
	nullify(r%next)
	return
	end subroutine newrecord


!!	Allocate a new tax record with default values

	subroutine newtax(t)
	type (taxlist), pointer :: t
	nullify(t)
	allocate(t)
	nullify(t%a)
	nullify(t%t)
	nullify(t%n)
	nullify(t%m)
	nullify(t%next)
	return
	end subroutine newtax

!!	Allocate a new netput with default values

	subroutine newnetput(np)
	type (netput), pointer :: np
	nullify(np)
	allocate(np)
	nullify(np%f)
	nullify(np%qref)
	nullify(np%pref)
	nullify(np%d)
	nullify(np%sd)
	nullify(np%t)
	np%nestcount = ' '
	np%nestfields = 0
	np%index = 0
	return
	end subroutine newnetput

!!	Allocate a new nest with default values

	subroutine newnest(n)
	type (nest), pointer :: n
	nullify(n)
	allocate(n)
	n%ident = ' '
	n%parentid = ' '
	n%value = ' '
	n%output = .false.
	nullify(n%f)
	nullify(n%next)
	nullify(n%parent)
        nullify(n%nlist)
	nullify(n%elements)
	nullify(n%subnests)
	n%index = 0
	n%level = 0
        n%assigned = .false.
	return
	end subroutine newnest

!!	Add d1 to union, excluding quoted singletons

	subroutine union(d1,u)
	type (domain), pointer :: d1, u, d, dd

	if (.not.associated(d1)) stop 'd1 is not associated'
	if (.not.associated(u)) nullify(u)

	d => d1 
	nullify(dd)
	if (associated(u)) dd => u
	do while (associated(d))
	  if (.not.singleton(d%s%name) .and. .not.indomain(d%s,u)) then
	    if (associated(u)) then
	      call newdomain(dd%next)
	      dd => dd%next
	    else
	      call newdomain(u)
	      nullify(dd)
	      dd => u
	    end if
	    dd%s => d%s
	  end if
	  d => d%next
	end do
	return
	end subroutine union

!!	Determine the subdomain of d1 in d2, i.e. those element of
!!	d1 which are not in d2:
!!>
!!                  subdomain = d1 \ (d1 ^ d2)
!!<
!!	The output is the set of elements in d1 which are not 
!!	singletons and not in d2.
!!
	subroutine findsubdomain(d1,d2,subdomain)
	type (domain), pointer :: d1, d2, subdomain, d, dd

	nullify(subdomain)

!	Look at each element in d1:

	d => d1 
	do while (associated(d))
	  if (.not.indomain(d%s,d2)) then
	    if (associated(subdomain)) then
	      call newdomain(dd%next)
	      dd => dd%next
	    else
	      call newdomain(subdomain)
	      dd => subdomain
	    end if
	    dd%s => d%s
	  end if
	  d => d%next
	end do
	return
	end subroutine findsubdomain

	logical function singleton(string)
	character (len=*) :: string
	singleton = ((string(1:1).eq.'"') .or. (string(1:1).eq."'"))
	return
	end function singleton

	subroutine writeparam(ident, index, d, string)
	character (len=*) :: ident, string
	type (domain), pointer :: d
	integer :: index
	character (len=255) :: domstr

	call writedomain(d, domstr)
	if (domstr.eq.' ') then
	  write(string,'(a,i3.3)') trim(ident), index
	else
	  write(string,'(a,i3.3,a)') trim(ident), index,'('//trim(domstr)//')'
	end if
	return
	end subroutine writeparam

	subroutine domainlabels(d,dl)
        type (domain), pointer :: d, dd	
	character (len=32) ::  dl(10)
	integer :: i
	dl = ' '
        dd => d
	i = 0
        do while (associated(dd))
	  i = i + 1
	  dl(i) = dd%s%name
          dd => dd%next
        end do
        return
        end subroutine domainlabels

	subroutine writedomain(d, string)
	character (len=*) :: string
	type (domain), pointer :: d, dd
	dd => d
	string = ' '
	do while (associated(dd))
	  if (string.eq.' ') then
	    string = dd%s%name
	  else
	    string = trim(string)//','//dd%s%name
	  end if
	  dd => dd%next
	end do
	return
	end subroutine writedomain

!!	Read through variable declarations

	subroutine readvar
	type (symbol), pointer :: s
	type (function), pointer :: fn
	type (field), pointer :: f
	type (record), pointer :: r
	character (len=9) :: vtype
	logical :: pfield
	integer :: fieldno

! Temporary KP
	write(*,*)'Calling Readvar'

	fn => functionlist
	do while (associated(fn))
	  r => fn%firstrecord
	  f => r%firstfield
	  write(*,fmt) trim(f%label)
          if ( f%label.eq.'$model') modelname = f%data
	  if (	f%label.ne.'$sectors'   .and. f%label.ne.'$commodities' .and. &
		f%label.ne.'$consumers' .and. f%label.ne.'$auxiliary') then
	    fn => fn%next
	    cycle 
	  else
	    select case(f%label)
	      case ('$sectors') 
	        vtype = 'sector'
	      case ('$commodities')
	        vtype = 'commodity'
	      case ('$consumers')
	        vtype = 'consumer'
	      case ('$auxiliary')
	        vtype = 'auxiliary'
	    end select
	  end if

	  if (log) write(*,fmt) ' ** '//trim(f%label)

	  if (r%nfield.ne.1)  stop 'Error: error in variable declaration definition record.'
	  f => r%firstfield
	  if (f%data.ne.' ') stop 'Error: trailing items not permitted with variable declaration title.'

	  do while (associated(r%next))
	    r => r%next
	    f => r%firstfield
	    if (r%nfield.ne.1)  stop 'Error: only one field per record in declarations.'
	    if (f%label.ne.' ') stop 'Error: labels are not expected in declaration recors.'
	    write(*,fmt)  vtype//' '//trim(f%data)
	    call findid(f%identtext, s)
	    if (associated(s)) stop 'Error: variable identifier has appeared previously.'
	    call addid(f%identtext, vtype, s, f%line)
	    f%ident => s

	    if(.not.associated(f%ident)) then
	      call InitializeMessages('all',0)
	      call InternalErrorStop('Readvar: f%ident is not associated.')
	    end if
!NO	    if(.not.associated(fn%firstrecord%firstfield%ident)) then
!NO	      call InitializeMessages('all')
!NO	      call InternalErrorStop('Readvar: fn%firstrecord%firstfield%ident is not associated.')
!NO	      write(*,*) 'fn%firstrecord%firstfield%ident added - KP'
!NO	      fn%firstrecord%firstfield%ident => s	      
!NO	    end if

	    s%dim = dimension(f%arg)
	    s%arg => f%arg   ! KP
	    call FindLabelFromLine(f%line, s%labelinfo) ! KP
	  end do
	  fn => fn%next
	end do
	return
	end subroutine readvar

!!	Read through records for each of the functions and add symbols
!!	to the master list as they are encountered.  Check that the 
!!	input data conforms to syntax rules.

	subroutine chkprod

	type (symbol), pointer :: s
	type (function), pointer :: fn
	type (field), pointer :: f
	type (record), pointer :: r
	logical :: pfield
	integer :: fieldno


	fn => functionlist

	do while (associated(fn))
	  r => fn%firstrecord
	  f => r%firstfield
	  if (f%label.ne.'$prod') then
	    fn => fn%next
	    cycle 
	  end if
	  if (log) write(*,fmt) ' ** '//trim(f%data)

!	Skip the elasticity nests on the first pass.

	  call InitializeMessages('pre',0)  ! KP
	  Messages%PreLine(1) = 'Processing function '//trim(f%label)//':'//trim(f%data) ! KP

	  do while (associated(r%next))

	    r => r%next
	    pfield = .false.	! Flag to assure only one P: field per record
	    f => r%firstfield
	    fieldno = 1

	    do while (associated(f))

	      if (log) write(*,fmt) ' -- '//trim(f%label)//':'//trim(f%data)

	      Messages%PreLine(2) = 'Processing field '//trim(f%label)//':'//trim(f%data) ! KP
	      call chkfield(f,fieldno,pfield)
!  Re-initialize 'pre' messages after number 1 after have finished this field ! KP
	      call InitializeMessages('pre',1) ! KP

	      fieldno = fieldno + 1
	      f => f%next
	    end do
	  end do
	  fn => fn%next
	end do

	call InitializeMessages('pre',0)  ! KP

	return
	end subroutine chkprod

!!	Evaluate fields in a production block record:
!!	- identifier dimension and type must be consistent
!!	- i: and o: must only appear as leading fields and reference prices
!!	- q: field must be second and reference a parameter 
!!	- p: fields may appear only once and must reference parameter
!!	- a: field must reference a consumer
!!	- n: field must reference an auxiliary variable 
!!	- m: field must reference a parameter
!!	- t: field must reference a parameter
!!
	subroutine chkfield(f,fieldno,pfield)
	integer :: fieldno
	logical :: pfield
	type (field), pointer :: f
	type (symbol), pointer :: s

!	We are scanning the fields assoiciated with named labels, hence we can skip
!	the nest-assignment fields:

	if (f%data.eq.' ' .or. f%data(1:1).eq.'$') then
	  if (index(f%label,'.tl').gt.0) stop ' *** Error: .tl syntax is not yet implemented.'
	  if (index(f%label,'.te').gt.0) stop ' *** Error: .te syntax is not supported.'
	  return
	endif

	select case (f%label)	! Evaluate fields which are not nesting assignments

	case('i')
	  if (fieldno.ne.1) stop ' *** Error: I: must be the leading field.'
	  call findid(f%identtext, s)
	  if (associated(s)) then 
	    if (s%type.ne.'commodity') stop ' *** Error: I: field must reference a commodity price'
	    if (s%dim.eq.-1) then
	      s%dim = dimension(f%arg)
	    else
	      if (s%dim.ne.dimension(f%arg))  stop ' *** Error: dimension of commodity price is different.'
	    end if
	    f%ident => s
	  else
! Change [30/03/04] by Ken. This is now an error since have processed all commodities in the $commodities part of the file
	    call ErrorStop('The price in this field is not declared in the $commodities section')
! NO 	    call addid(f%identtext, 'commodity', s, f%line)
! NO	    f%ident => s
! NO	    s%dim = dimension(f%arg)
	  end if

	case('o')
! Tom	  if (fieldno.ne.1) stop ' *** Error: O: must be the leading field.'
! KP change to use ErrorStop routine
	  if (fieldno.ne.1) call ErrorStop('O: must be the leading field.')
	  call findid(f%identtext, s)
	  if (associated(s)) then 
! (was. KP) if (s%type.ne.'commodity') stop ' *** Error: O: field must reference a commodity price'
	    if (s%type.ne.'commodity') call ErrorStop('O: field must reference a commodity price')
	    if (s%dim.eq.-1) then
	      s%dim = dimension(f%arg)
	    else
! KP	      if (s%dim.ne.dimension(f%arg))  stop ' *** Error: dimension of commodity price is different.'
	      if (s%dim.ne.dimension(f%arg))  call ErrorStop('dimension of commodity price is different from its'// &
			'\n declaration in the $commodities section.')
	    end if
	    f%ident => s
	  else
! Change [30/03/04] by Ken. This is now an error since have processed all commodities in the $commodities part of the file
	    call ErrorStop('The price in this field is not declared in the $commodities section')
! NO	    call addid(f%identtext, 'commodity', s, f%line)
! NO	    f%ident => s
! NO	    s%dim = dimension(f%arg)
	  end if

	case('q')
	  if (fieldno.ne.2) stop ' *** Error: Q: must be in the second field.'
	  if (f%identtext(1:1).eq.'(') return

! [31/03/04] Allow "-" followed by parameter in q: field ! KP
	  if (f%identtext(1:1).eq.'-') return ! KP

	  if (numeric(f%data)) return
	  call findid(f%identtext, s)
	  if (associated(s)) then 
	    if (s%type.ne.'parameter') stop ' *** Error: Q: field must reference a parameter'
	    if (s%dim.eq.-1) then
	      s%dim = dimension(f%arg)
	    else
	      if (s%dim.ne.dimension(f%arg))  stop ' *** Error: dimension of parameter is different.'
	    end if
	    f%ident => s
	  else
	    call addid(f%identtext, 'parameter', s, f%line)
	    f%ident => s
	    s%dim = dimension(f%arg)
	  end if

	case('p')
	  if (pfield) stop ' *** Error: P: field appears more than once.'
	  pfield = .true.
	  if (f%identtext(1:1).eq.'(') return
	  if (numeric(f%data)) return
	  call findid(f%identtext, s)
	  if (associated(s)) then 
	    if (s%type.ne.'parameter') stop ' *** Error: P: field must reference a parameter'
	    if (s%dim.eq.-1) then
	      s%dim = dimension(f%arg)
	    else
	      if (s%dim.ne.dimension(f%arg))  stop ' *** Error: dimension of parameter is different.'
	    end if
	    f%ident => s
	  else
	    call addid(f%identtext, 'parameter', s, f%line)
	    f%ident => s
	    s%dim = dimension(f%arg)
	  end if

	case('a')
	  call findid(f%identtext,s)
	  if (associated(s)) then 
	    if (s%type.ne.'consumer') stop ' *** Error: A: field must reference a consumer'
	    if (s%dim.eq.-1) then
	      s%dim = dimension(f%arg)
	    else
	      if (s%dim.ne.dimension(f%arg))  stop ' *** Error: dimension of parameter is different.'
	    end if
	    f%ident => s
	  else
	    call addid(f%identtext, 'consumer', s, f%line)
	    f%ident => s
	    s%dim = dimension(f%arg)
	    s%arg => f%arg   ! KP
	    call FindLabelFromLine(f%line, s%labelinfo) ! KP
	  end if


	case('n')
	  call findid(f%identtext,s)
	  if (associated(s)) then 
	    if (s%type.ne.'auxiliary') stop ' *** Error: N: field must reference an AUXILIARY'
	    if (s%dim.eq.-1) then
	      s%dim = dimension(f%arg)
	    else
	      if (s%dim.ne.dimension(f%arg))  stop ' *** Error: dimension of AUXILIARY variable changes.'
	    end if
	    f%ident => s
	  else
	    call addid(f%identtext, 'auxiliary', s, f%line)
	    f%ident => s
	    s%dim = dimension(f%arg)
	  end if

	case('m')
	  if (f%identtext(1:1).eq.'(') return
	  if (numeric(f%data)) return
	  call findid(f%identtext,s)
	  if (associated(s)) then 
	  if (s%type.ne.'parameter') stop ' *** Error: M: multiplier must reference an parameter'
	    if (s%dim.eq.-1) then
	      s%dim = dimension(f%arg)
	    else
	      if (s%dim.ne.dimension(f%arg))  stop ' *** Error: dimension of parameter variable changes.'
	    end if
	    f%ident => s
	  else
	    call addid(f%identtext, 'parameter', s, f%line)
	    f%ident => s
	    s%dim = dimension(f%arg)
	  end if

	case('t')
	  if (f%identtext(1:1).eq.'(') return
	  if (numeric(f%data)) return
	  call findid(f%identtext,s)
	  if (associated(s)) then 
	    if (s%type.ne.'parameter') stop ' *** Error: T: multiplier must reference an parameter'
	    if (s%dim.eq.-1) then
	      s%dim = dimension(f%arg)
	    else
	      if (s%dim.ne.dimension(f%arg))  stop ' *** Error: dimension of parameter variable changes.'
	    end if
	    f%ident => s
	  else
	    call addid(f%identtext, 'parameter', s, f%line)
	    f%ident => s
	    s%dim = dimension(f%arg)
	  end if

	case default
	  stop ' *** Error: unrecognized field label in $prod block'

	end select 

	return

	end subroutine chkfield

!!	Read nest assignments from the first record of each $prod block

	subroutine readnest

	type(function), pointer :: fn
	type(record), pointer :: r
	type(field), pointer :: f
	type(nest), pointer ::	esub, &	! Top level input nest
				etrn, &	! Top level output nest
				n, &	! Last nest allocated
				nn,np	! Nest pointer

	logical :: subnest	! Logical flag that a subnest has been specified.
	integer :: fieldno
	character (len=255) :: nestid

!	First pass to set up nesting pointers:

	fn => functionlist
	nullify(f)
	nullify(n)
	do while (associated(fn))
	  r => fn%firstrecord
	  f => r%firstfield

!	Set up the top level input and output nests:

	  call newnest(esub)
	  fn%inputs  => esub
	  esub%value = '0'
	  esub%ident = 's'

	  call newnest(etrn)
	  fn%outputs => etrn
	  etrn%value = '0'
	  etrn%ident = 't'
	  etrn%output = .true.

	  esub%next => etrn
	  n => etrn
	  subnest = .false.
	  fieldno = 1

!	Look through all the fields on this record following the first:

	  do while (associated(f%next))
	    f => f%next
	    fieldno = fieldno + 1
	    if (f%spantext.ne.' ') stop ' *** Span notation "#" not permited in elasticity.'

!	If an elasticity is specified as a parameter, verify that the parameter domain
!	is a subset of the function domain:

	    if ((f%argtext.ne.' ') .and. (.not.subset(f%arg, r%firstfield%arg))) then
    	      write(*,fmt) ' '//trim(r%firstfield%label)//':'//trim(r%firstfield%data)//' => '//trim(r%firstfield%argtext)
    	      write(*,fmt) ' '//trim(f%label)//':'//trim(f%data)//' => '//trim(f%argtext)
	      stop ' ***  Invalid elasticity field domain.'
	    end if

!	Nest ID is specfied in the field label:

    	    nestid = f%label
	    if (index(nestid,'.tl').gt.0) stop ' *** Error: .tl syntax is not yet implemented.'
	    if (index(nestid,'.te').gt.0) stop ' *** Error: .te syntax is not supported.'

	    if (nestid.eq.'s') then
	      if (subnest) stop ' *** s: field must preceed subnest fields.'
	      if (associated(esub%f)) stop ' *** Multiple s: fields in $prod block'
	      if (fieldno.gt.3)       stop ' *** s: field must be first or second elasticity'
	      esub%f => f
	      esub%value = f%data
	      cycle
	    endif

	    if (nestid.eq.'t') then
	      if (subnest) stop ' *** t: field must preceed subnest fields.'
	      if (associated(etrn%f)) stop ' *** Multiple t: fields in $prod block'
	      if (fieldno.gt.3)       stop ' *** t: field must be first or second elasticity'
	      etrn%f => f
	      etrn%value = f%data
	      cycle
	    endif

	    subnest = .true.
	    call newnest(n%next)
	    n => n%next

	    call parsenestid(nestid,n%ident,n%parentid)
	    n%f => f
	    n%value = f%data
	  end do

!	Now install pointers to parents where an explicit parent has
!	been declared:

	  n => fn%inputs
	  do while (associated(n))
	    if (n%parentid.eq.' ') then
	      if (.not.(associated(n,fn%inputs) .or. associated(n,fn%outputs))) n%level = 1
	      n => n%next
	      cycle
	    end if

!	An explicit parent has been declared - find it

	    call findnest(n%parentid, fn, np)
	    n%parent => np
	    n%level = np%level + 1

!	Add this nest to the list of subnests in p:

	    call addsubnest(n, np)
	    n => n%next
	  end do 
	    
          if (log .and. associated(fn%inputs)) then
!.	    write(*,fmt) ' Nests:'
	    n => fn%inputs
	    do while (associated(n))
!.	      if (associated(n%parent)) write(*,fmt) ' '//trim(n%ident)//' -> '//trim(n%parentid)
	      n => n%next
	    end do
	  end if

	  fn => fn%next
	end do

	return
	end subroutine readnest

!!	Split a nest elasticity field into a nest identifier and a 
!!	parent identifier:

	subroutine parsenestid(nestid,ident,parentid)
	character (len=*) :: nestid	!! Identifier for nest
	character (len=*) :: ident	!! Identifier of the nest itself
	character (len=*) :: parentid	!! Identifier of the nest parent 

	integer :: ic, jc
	ic = index(nestid,'(')
	jc = index(nestid,')',back=.true.)
	if (ic.gt.0 .and. jc.eq.0) then
	  write(*,fmt) trim(nestid)
	  stop ' Invalid parent node -- missing closing parenthesis.'
	end if
	if (ic.eq.0 .and. jc.gt.0) then
	  write(*,fmt) trim(nestid)
	  stop ' Invalid parent node -- missing opening parenthesis.'
	end if
	if (ic.eq.0) then
	  ident = nestid
	  parentid = ' '
	else
	  if (jc.lt.ic+2) then
	    write(*,fmt) trim(nestid)
	    stop ' Invalid parent node -- parenthesis mismatch.'
	  end if
	  ident = nestid(1:ic-1)
	  parentid = nestid(ic+1:jc-1)
	end if
	return
	end subroutine parsenestid

!!	Read through the production function coefficients:
!!	- Create netput descriptions of each production input and output
!!	- Create linked list data structure characterizing the nested CES function

	subroutine readcoef
	type (function), pointer :: fn
	type (symbol), pointer :: s
	type (record), pointer :: r
	type (netput), pointer :: np
	type (domain), pointer :: d
	type (field), pointer :: f, fnf, taxagent, lastfield
	type (nest), pointer :: n
	logical :: input, conditional
	integer :: nparent

	nullify(taxagent)	! No tax agent defined at this point
	  
	fn => functionlist
	do while (associated(fn))
	  r => fn%firstrecord
	  fnf => r%firstfield
	  if (fnf%label.ne.'$prod') then
	    fn => fn%next
	    cycle
	  end if

! Extract field identifier and add it to the field record - added by KP (?omitted by TR?)
	  call findid(fnf%identtext, s)

! 29/01/04 This is probably unnecessary becase of addition of f%ident in routine parse
! Add this identifier s to the field record - added by KP (?omitted by TR?)
	  if(.not.(associated(fnf%ident))) then
	    write(*,*)'KP add - was not already asociated'
	    fnf%ident => s
	  else
	    write(*,*)'KP add - was already associated'
	  end if

!	Read all the records in this function after the elasticity record:

	  do while (associated(r%next))
	    r => r%next

!	Look through the fields in this record:

	    f => r%firstfield
	    nparent = 0
	    call findid(f%identtext, s)

!	Create a netput for this record:

	    call newnetput(np)
	    np%f => f
	    np%output = (f%label.eq.'o')

	    np%nestcount = '0'
	    np%nestfields = 0
	    conditional = .true.

	    if (associated(fnf%arg)) call union(fnf%arg, np%d)
	    if (associated(f%arg))   call union(f%arg,  np%d)
	    if (associated(f%span))  call union(f%span, np%d)
	    if (associated(np%d))    call findsubdomain(np%d, fnf%arg, np%sd)
	    nullify(lastfield)	! Reading the first field -- no prior field at this point


!	Look for nest identifiers:

	    do while (associated(f%next))
	      f => f%next

!	A nest identifier is a field with no data component:

	      if (f%data.eq.' ') then

!	Locate this nest identifier in the master list:

	        call findnest(f%label, fn, n)

!	Record the nest type:

		call nesttype(fn, n, np%output)

!	Add this netput to this nest:

		call addnetput(np, n, f)
		nullify(lastfield)
		np%nestfields = np%nestfields + 1
		if (f%condtn.ne.' ') then
		  np%nestcount = trim(np%nestcount)//'+1$'//trim(f%condtn)
		else
		  np%nestcount = trim(np%nestcount)//'+1'
		  conditional = .false.
		end if
	      else
	        call readelem(np, f, taxagent, lastfield)
		lastfield => f
	      end if
	    end do

!	Assign this netput to a toplevel

	    if (     np%output) n => fn%outputs
	    if (.not.np%output) n => fn%inputs
	    if (np%nestfields.eq.0) then
	      nullify(f)
	      call addnetput(np, n, f)
	    else
	      if (conditional) then
	        call newfield(f)
		f%condtn = '('//trim(np%nestcount)//'=0)'
		call addnetput(np, n, f)
	      endif
	    end if
	  end do
	  fn => fn%next
	end do

	return
	end subroutine readcoef

!!	Introduce type definitions for nests beginning at n
!!	and proceeding to the top of the tree.

	subroutine nesttype(fn, n, output)
	type (function), pointer :: fn	!! Function we fare looking at
	type (nest), pointer :: n	!! Nest from which to trace
	logical :: output		!! Flag for whether this is an output nest
	type (nest), pointer :: np
	if (n%output .and. .not.output) stop 'Error: input enters an output nest'
	n%output = output
	if (.not. associated(n%parent)) then
	  if (output) then
	    n%parent => fn%outputs
	    call addsubnest(n,fn%outputs)
	  else
	    n%parent => fn%inputs
	    call addsubnest(n,fn%inputs)
	  end if
	  return
	else
	  np => n
	  do while (associated(np%parent))
	    np => np%parent
	    if (np%output .and. .not.output) stop 'Error: input enters an output nest'
	    np%output = output
	  end do
	  if (output      .and. .not.associated(np,fn%outputs)) then
	    np%parent => fn%outputs
	    call addsubnest(np%parent,fn%outputs)
	  end if
	  if (.not.output .and. .not.associated(np,fn%inputs))  then
	    np%parent => fn%inputs
	    call addsubnest(np%parent,fn%inputs)
	  end if
	end if
	return
	end subroutine nesttype

!!	Read a single field and introduce the associated logic
!!	in a netput descriptor:

	subroutine readelem(np, f, taxagent, lastfield)
	type (netput), pointer :: np	!! Netput to which this should be assigned
	type (field), pointer  :: f	!! Field for which logic is to be introduced
	type (field), pointer :: taxagent	!! Field which defines applicable tax agent field
	type (field), pointer :: lastfield	!! Pointer to the previous field in this record (if one)
	type (taxlist), pointer :: t

	select case (f%label)

	case ('p')
	  if (associated(np%pref)) stop 'P: field appears more than once.'
	  np%pref => f

	case ('q')
	  if (associated(np%qref)) stop 'Q: field appears more than once.'
	  np%qref => f

	case ('a')
	  taxagent => f

	case ('t')
	  if (.not.associated(taxagent)) stop 'Tax field appears prior to specification of a tax agent (A:).'
	  call addtax(np,t)
	  t%t => f
	  t%a => taxagent

	case ('n')
	  if (.not.associated(taxagent)) stop 'Tax field appears prior to specification of a tax agent (A:).'
	  call addtax(np,t)
	  t%n => f
	  t%a => taxagent
	  
	case ('m')
	  if (.not.associated(taxagent)) stop 'Tax field appears prior to specification of a tax agent (A:).'
	  call lasttax(np,t)
	  if (.not.associated(lastfield))     stop 'M: field must follow N: field'
	  if (.not.associated(t%n,lastfield)) stop 'M: field must follow immediately after N: field'
	  t%m => f

	end select

	lastfield => f
	return
	end subroutine readelem

!!	Create symbol table entries for the production sectors

	subroutine sectors

	character (len=31) :: label
	type(function), pointer :: fn
	type(record), pointer :: r
	type(field), pointer :: f
	type(symbol), pointer :: s
	type(domain), pointer :: d

        fn => functionlist
        nullify(f);
        do while (associated(fn))
	  r => fn%firstrecord
	  f => r%firstfield
	  if (log) write(*,fmt) ' '//trim(f%label)//':'//trim(f%data)
	  fn => fn%next
        end do

	nullify(f);
	fn => functionlist
	do while (associated(fn))
	  r => fn%firstrecord
	  f => r%firstfield
	  if (log) write(*,fmt) ' label: '//trim(f%label)

!  Re-initialize 'pre' messages again at start of each function  ! KP
	  call InitializeMessages('pre',0)
	  Messages%PreLine(1) = 'Processing function '//trim(f%label)//':'//trim(f%data)

	  if (f%label.eq.'$prod') then
	    if (log) write(*,fmt) ' Checking $prod for: '//trim(f%identtext) !KP
	    call findid(f%identtext, s)
	    if (.not.associated(s)) then
	      call ErrorStop('The activity level '//trim(f%identtext)//' in this $prod block is not declared in the $sectors section')
	    end if
! KP addition ??
	    if(.not.associated(f%ident)) then
		write(*,fmt) '$prod. f%identtext=', trim(f%identtext), ' not associated - added by KP'
	        f%ident => s
	    end if

	    if (s%dim.ne.dimension(f%arg)) call ErrorStop('$prod variable dimension is inconsistent')
	    if (f%spantext.ne.' ') call ErrorStop('Span notation "#" not permited in $prod block.')
	    d => f%arg
	    do while (associated(d))
	      if (d%s%name(1:1).eq.'"' .or. d%s%name(1:1).eq."'") &
		call ErrorStop('Singleton sets not permitted in $prod domain')
	      d => d%next
	    end do
	  end if

          if (f%label.eq.'$demand') then
            if (log) write(*,fmt) ' Checking $demand for '//trim(f%identtext) ! KP
            call findid(f%identtext, s)
            if (.not.associated(s)) then
	      call ErrorStop('The revenue '//trim(f%identtext)//' in this $demand block is not declared in the $consumers section')
            end if

! KP addition ??
	    if(.not.associated(f%ident)) then
		write(*,fmt) '$demand. f%identtext=', trim(f%identtext), ' not associated - added by KP'
	        f%ident => s
	    end if

	    if(log) write(*,fmt)'Sectors. $demand. f%identtext is:', trim(f%identtext),' FINDID ok' ! KP
            if (s%dim.ne.dimension(f%arg)) call ErrorStop('$demand variable dimension is inconsistent')
            if (f%spantext.ne.' ') call ErrorStop('Span notation "#" not permited in $demand block.')
            d => f%arg
            do while (associated(d))
              if (d%s%name(1:1).eq.'"' .or. d%s%name(1:1).eq."'") &
                call ErrorStop('Singleton sets not permitted in $demand domain')
              d => d%next
            end do
          end if

! Addition [02/04/04] by KP to check variable in $constraint: line
          if (f%label.eq.'$constraint') then
            if (log) write(*,fmt) ' Checking $constraint for '//trim(f%identtext)
            call findid(f%identtext, s)
            if (.not.associated(s)) then
	      call ErrorStop('The variable '//trim(f%identtext)//' in this $constraint line is not declared in the $auxiliary section')
            end if

! KP addition
	    if(.not.associated(f%ident)) then
		write(*,fmt) '$constraint. f%identtext=', trim(f%identtext), ' not associated - added by KP'
	        f%ident => s
	    end if

            if (s%dim.ne.dimension(f%arg)) call ErrorStop('$constraint variable dimension is inconsistent')
            if (f%spantext.ne.' ') call ErrorStop('Span notation "#" not permited in $constraint block.')
            d => f%arg
            do while (associated(d))
              if (d%s%name(1:1).eq.'"' .or. d%s%name(1:1).eq."'") &
                call ErrorStop('Singleton sets not permitted in $constraint domain') ! Is this true? (KP)
              d => d%next
            end do
          end if

	  fn => fn%next
	end do
	return
	end subroutine sectors

!!	Determine which a string is numeric:

	logical function numeric(string)
	character (len=*) :: string	!! String to be evaluated (.true. if a Fortran-readable numeric value)
	real :: testval
	integer :: stat
	if (index('.+-0123456789', string(1:1)).gt.0) then
	  read(string,*,iostat=stat) testval
!	  if (stat.ne.0) stop ' ** Error in numeric field.'
	  if (stat.ne.0) call ErrorStop( ' ** Error in numeric field.')
	  numeric = .true.
	else
	  numeric = .false.
	end if
	return
	end function numeric

!!	Add a netput to a nest.

	subroutine addnetput(np, n, f)
	type (netput), pointer :: np	!! Netput to be added
	type (nest), pointer :: n	!! Nest in which the netput is to be introduced
	type (field), pointer :: f	!! Field in which the nest assignment is made (for condition if any)
	type (netputlist), pointer :: nl
	if (.not.associated(n%elements)) then
	  nullify(nl)
	  allocate(nl)
	  n%elements => nl
	else
	  nl => n%elements
	  do while (associated(nl%next))
	    nl => nl%next
	  end do
	  allocate(nl%next)
	  nl => nl%next
	endif
	nullify(nl%np)
	nullify(nl%next)
	nl%np => np
	if (associated(f)) nl%f => f
	return
	end subroutine addnetput

!!	Add a tax coefficient to the associated linked list for a production netput:

	subroutine addtax(np,t)
	type (netput), pointer :: np	!! Netput to which tax is applied
	type (taxlist), pointer :: t	!! Pointer to the tax field

	if (associated(np%t)) then
	  call lasttax(np, t)
	  call newtax(t%next)
	  t => t%next
	else
	  call newtax(np%t)
	  t => np%t
	end if
	return
	end subroutine addtax

!!	Set pointer to the last tax applied to a given netput:

	subroutine lasttax(np, t)
	type (netput), pointer :: np	!! Netput to be queried
	type (taxlist), pointer :: t	!! On output points to the last tax coefficient

	t => np%t
	do while (associated(t%next))
	  t => t%next
	end do
	return
	end subroutine lasttax


!!	Generate a pointer to the nest with a given name in a particular function:

	subroutine findnest(nestid, fn, n)
	character (len=*) :: nestid	!! String descriptor of the nest to be located
	type (function), pointer :: fn	!! Function in which to locate the nest
	type (nest), pointer :: n	!! Pointer to the located nest

	n => fn%inputs
	if (log) write(*,fmt) ' Locating nest: '//trim(nestid)
	do while (associated(n))
	  if (nestid.eq.n%ident) return
	  n => n%next
	end do
	stop 'nestid not found'
	end subroutine findnest
	    
!!	At a subnest to a nest.  This requires creation of a nestlist if
!!	one does not already exist

	subroutine addsubnest(n, np)
	type (nest), pointer :: n	!! Nest to be added as a subnest
	type (nest), pointer :: np	!! Nest to which the subnest is to be introduced
	type (nestlist), pointer :: nl

	if (.not.associated(np%subnests)) then
	  nullify(nl)
	  allocate(nl)
	  np%subnests => nl
	  nullify(nl%next)
	  nl%n => n
          n%nlist => nl
	else
	  nl => np%subnests
	  do while (associated(nl%next))
	    nl => nl%next
	  end do
	  allocate(nl%next)
	  nl => nl%next
	  nullify(nl%next)
	  nl%n => n
          n%nlist => nl
	end if
	return
	end subroutine addsubnest

!!	Returns .true. if the argument is quoted text, i.e. a character
!!	string beginning and ending with a pair of single or double quotes.

	logical function quoted(arg)
	character (len=*) :: arg
	integer :: lc
	lc = len_trim(arg)
	quoted = .true.
	if (arg(1:1).eq."'" .and. arg(lc:lc).eq."'") return
	if (arg(1:1).eq.'"' .and. arg(lc:lc).eq.'"') return
	quoted = .false.
	return
	end function quoted

!!	Public function determines whether two strings are
!!	identical apart from case.  

        logical function eqv(t1,t2)
	implicit none
        integer :: i
        character (len=*) :: t1,t2
	character (len=255) :: s1,s2
	character (len=26) :: upper = 'ABCDEFGHIJKLMNOPQRSTUVWXYZ', &
			      lower = 'abcdefghijklmnopqrstuvwxyz'

	s1 = t1
	s2 = t2
	if (quoted(t1)) s1 = t1(2:len_trim(t1)-1)
	if (quoted(t2)) s2 = t2(2:len_trim(t2)-1)
	  
	eqv = .false.
	if (len_trim(s1).ne.len_trim(s2)) return
        do i=1,len_trim(s1)
          if (.not. ((s1(i:i).eq.s2(i:i)) .or. &
	      (index(upper,s1(i:i))+index(upper,s2(i:i)).ne.0 .and. &
	       index(upper,s1(i:i)).eq.index(lower,s2(i:i)) .and. &
	       index(upper,s2(i:i)).eq.index(lower,s1(i:i))))) return
	end do
	eqv = .true.
	return
	end function eqv

        subroutine consumers

        character (len=31) :: label
        type(function), pointer :: fn
        type(record), pointer :: r
        type(field), pointer :: f
        type(symbol), pointer :: s
        type(domain), pointer :: d

        nullify(f);
        fn => functionlist
        do while (associated(fn))
          r => fn%firstrecord
          f => r%firstfield
          fn => fn%next
        end do
        return
        end subroutine consumers


!!      Read through records for each of the demand functions and add symbols
!!      to the master list as they are encountered.  Check that the 
!!      input data conforms to syntax rules.

        subroutine chkdemand

        type (symbol), pointer :: s
        type (function), pointer :: fn
        type (field), pointer :: f,fnf
        type (record), pointer :: r
	type (endowmentlist), pointer :: el
	type (endowment), pointer :: e
        integer :: fieldno
	type (domain), pointer   :: d           ! Pointer to the domain of a D: field


        fn => functionlist
        do while (associated(fn))
          r => fn%firstrecord
          fnf => r%firstfield
          if (fnf%label.ne.'$demand') then
            fn => fn%next
            cycle 
          end if
	  f => r%firstfield
          if (log) write(*,fmt) ' ** '//trim(f%data)

	  call InitializeMessages('pre',0)  ! KP
	  Messages%PreLine(1) = 'Processing function '//trim(fnf%label)//':'//trim(fnf%data) ! KP
	  if(log) write(*,*) trim(Messages%PreLine(1))

          do while (associated(r%next))
            r => r%next
            f => r%firstfield

	    Messages%PreLine(2) = 'Processing field '//trim(f%label)//':'//trim(f%data) ! KP
	    if(log) write(*,*) trim(Messages%PreLine(2))

            fieldno = 1
	    if (f%label.eq.'d') then
	      if (associated(fn%demand)) stop '*** Error: multiple D: fields encountered.'
	      fn%demand => f
	      if (associated(f%arg)) then
	        if (associated(fnf%arg)) then
	          if (.not.subset(f%arg,fnf%arg)) stop '*** Error: D: fields must refer to single commodities.'
		else
		  nullify(d)
		  call union(f%arg,d)
		  if (d%dim.gt.0) stop '*** Error: D: field must refer to single commodities.'
		end if
	      end if
	      if (r%nfield.gt.2)  stop '*** Error: D: fields must be singleton fields, one per consumer.'
	    end if

	    if (f%label.eq.'e') then
	      if (.not.associated(fn%endowments)) then
	        allocate(el)	      
		fn%endowments => el
	      else
	        allocate(el%next)
		el => el%next
	      endif
	      nullify(el%next)
	      allocate(e)
	      nullify(e%q);nullify(e%r);nullify(e%d);nullify(e%sd)
	      el%e => e
	      e%f => f
	      if (associated(fnf%arg)) call union(fnf%arg, e%d)
	      if (associated(f%arg))   call union(f%arg,   e%d)
	      if (associated(f%span))  call union(f%span,  e%d)
	      if (associated(e%d))     call findsubdomain(e%d, fnf%arg, e%sd)
	    else
	      nullify(e)
	    end if

! Addition by KP to check that price in d: or e: field is declared in $commodities section
	    if (f%label.eq.'e' .or. f%label.eq.'d') then
	      call findid(f%identtext, s)
	      if (associated(s)) then 
	        if (s%type.ne.'commodity') then
		  call ErrorStop(f%label//' field must reference a price declared in the $commodities section')
	        end if
	        if (s%dim.eq.-1) then
	          s%dim = dimension(f%arg)
	        else
	          if (s%dim.ne.dimension(f%arg))  then
		    call ErrorStop( &
		      'Dimension of price in this field is different'// &
		      '\n from the dimension as declared in the $commodities section')
		  end if
	        end if
	        f%ident => s
	      else
! This is now an error since have processed all commodities in the $commodities part of the file
	        call ErrorStop('The price in this field is not declared in the $commodities section')
	      end if
	    end if ! if (f%label.eq.'e' .or. f%label.eq.'d') then

            do while (associated(f))
              if (log) write(*,fmt) ' -- '//trim(f%label)//':'//trim(f%data)
	      Messages%PreLine(2) = 'Processing field '//trim(f%label)//':'//trim(f%data) ! KP
              call chkdfield(f,fieldno,e)
!  Re-initialize 'pre' messages after number 1 after have finished this field ! KP
	      call InitializeMessages('pre',1) ! KP
              fieldno = fieldno + 1
              f => f%next
            end do
          end do
	  if (.not.associated(fn%demand)) stop '*** Error: require exactly one D: field per consumer.'
          fn => fn%next
        end do
        return
        end subroutine chkdemand

!!      Evaluate fields in a demand block record:
!!      - identifier dimension and type must be consistent
!!      - e: and d: must only appear as leading fields 
!!      - q: field must be second and reference a parameter 
!!      - r: field must reference an auxiliary and be the third
!!
        subroutine chkdfield(f,fieldno,e)
        integer :: fieldno
        type (field), pointer :: f
        type (symbol), pointer :: s
	type (endowment), pointer :: e

        select case (f%label)   ! Evaluate fields which are not nesting assignments

        case('d')
          if (fieldno.ne.1) stop ' *** Error: D: must be the leading field.'
          call findid(f%identtext, s)
          if (associated(s)) then 
            if (s%type.ne.'commodity') stop ' *** Error: D: field must reference a commodity price'
            if (s%dim.eq.-1) then
              s%dim = dimension(f%arg)
            else
              if (s%dim.ne.dimension(f%arg))  stop ' *** Error: dimension of commodity price is different.'
            end if
            f%ident => s
          else
	    call ErrorStop('The price in this d: field is not declared in the $commodities section')
          end if

        case('e')
          if (fieldno.ne.1) stop ' *** Error: E: must be the leading field.'
          call findid(f%identtext, s)
          if (associated(s)) then 
            if (s%type.ne.'commodity') stop ' *** Error: E: field must reference a commodity price'
            if (s%dim.eq.-1) then
              s%dim = dimension(f%arg)
            else
              if (s%dim.ne.dimension(f%arg))  stop ' *** Error: dimension of commodity price is different.'
            end if
            f%ident => s
          else
! [02/04/04] Change by KP. This is now an error since have processed all commodities in the $commodities part of the file.
	    call ErrorStop('The price in this e: field is not declared in the $commodities section')
          end if

        case('q')
	  if (.not.associated(e)) stop '*** Error: Q: field only expected in an endowment record with a $demand block'
          if (fieldno.ne.2) stop ' *** Error: Q: must be in the second field.'
	  e%q => f
          if (f%identtext(1:1).eq.'(') return

! [31/03/04] Allow "-" followed by parameter in q: field ! KP
	  if (f%identtext(1:1).eq.'-') return ! KP

          if (numeric(f%data)) return
          call findid(f%identtext, s)
          if (associated(s)) then 
            if (s%type.ne.'parameter') stop ' *** Error: Q: field must reference a parameter'
            if (s%dim.eq.-1) then
              s%dim = dimension(f%arg)
            else
              if (s%dim.ne.dimension(f%arg))  stop ' *** Error: dimension of parameter is different.'
            end if
            f%ident => s
          else
            call addid(f%identtext, 'parameter', s, f%line)
            f%ident => s
            s%dim = dimension(f%arg)
          end if

        case('r')
	  if (.not.associated(e)) stop '*** Error: R: field only expected in an endowment record with a $demand block'
	  if (associated(e%r)) stop '*** Error: multiple R: fields?'
	  e%r => f
          if (f%identtext(1:1).eq.'(')	stop ' *** Error: R: field must reference an auxiliary'
          if (numeric(f%data))		stop ' *** Error: R: field must reference an auxiliary'
          call findid(f%identtext, s)
          if (associated(s)) then 
            if (s%type.ne.'auxiliary') stop ' *** Error: R: field must reference an auxiliary'
            if (s%dim.eq.-1) then
              s%dim = dimension(f%arg)
            else
              if (s%dim.ne.dimension(f%arg))  stop ' *** Error: dimension of auxiliary is different.'
            end if
            f%ident => s
          else
	    call ErrorStop('The variable in this r: field is not declared in the $auxiliary section')
          end if

        case default
          stop ' *** Error: unrecognized field label in $prod block'

        end select 

        return

        end subroutine chkdfield

	subroutine UserOrInternalErrorStop(ErrorType, message)  ! KP
! Stops with fatal error. Give message
! Input
	character (len=*) ErrorType ! 'u' (user) or 'i' (internal)
	character (len=*) message
! Local
	character (len=255) :: OrigPreLine1, MessageLine
	character (len=2048) :: MessageLeft
	integer :: n1
	logical :: MessageAllWritten

! Store Messages%PreLine(1)
	OrigPreLine1=Messages%PreLine(1)
! Modify Messages%PreLine(1) - but only if it is not all blank
	if(OrigPreLine1.ne.' ') then
	  Messages%PreLine(1)='Error while '//trim(OrigPreLine1)
	else
	  Messages%PreLine(1)=' '
	end if
! Put out the relevant PreLines
	write(*,*) ' '
	call PrePostErrorWarning(Size(Messages%PreLine), Messages%PreLine)
	if(ErrorType.eq.'u') then
	  write(*,fmt)'Fatal error. Message follows.'
	else
	  write(*,fmt)'Internal program error. Message follows.'
	end if

! Output the message
! "\n" in Message means start new line
	MessageLeft = Message
	MessageAllWritten = .false.
	do while(.not.MessageAllWritten)
	  n1 = INDEX(MessageLeft, '\n')
	  if(n1.gt.0) then
	    write(*,'(2x,a)') trim(MessageLeft(1:n1-1))
	    MessageLeft = MessageLeft(n1+2:)
	  else
	    write(*,'(2x,a)') trim(MessageLeft)
	    MessageAllWritten = .true.
	  end if
	end do ! while(.not.MessageAllWritten)

! Put out the relevant PostLines
	call PrePostErrorWarning(Size(Messages%PostLine), Messages%PostLine)
	if(ErrorType.eq.'i') then
	  write(*,fmt)'** Please tell the developers of this code. **'
	end if
	write(*,*) ' '
!
!	stop
!  Next sets DOS errorlevel 1 and stops. This works with LF95
	call EXIT(1)
!
	end subroutine UserOrInternalErrorStop

	subroutine ErrorStop(message)  ! KP
! This is for user errors (as distinct from internal program errors)
! Stops with fatal error. Give message
! Input
	character (len=*) message
!
	call UserOrInternalErrorStop('u', message)
!
	return
	end subroutine ErrorStop

	subroutine InternalErrorStop(message)  ! KP
! This is for internal program errors (as distinct from user errors)
! Stops with fatal error. Give message
! Input
	character (len=*) message
!
	call UserOrInternalErrorStop('i', message)
!
	return
	end subroutine InternalErrorStop

	subroutine GiveWarning(message)  ! KP
! Gives a warning.
! Input
	character (len=*) message
! Local
	character (len=255) :: OrigPreLine1

! Store Messages%PreLine(1) so can restore it after modifying it
	OrigPreLine1=Messages%PreLine(1)
! Modify Messages%PreLine(1)
	Messages%PreLine(1)='Warning while '//trim(OrigPreLine1)
! Put out the relevant PreLines
	write(*,*) ' '
	call PrePostErrorWarning(Size(Messages%PreLine), Messages%PreLine)
	write(*,fmt) ' %% Warning. Message follows.'
	write(*,fmt) message
! Put out the relevant PostLines
	call PrePostErrorWarning(Size(Messages%PostLine), Messages%PostLine)
	write(*,*) ' '
! Restore original PreLine(1)
	Messages%PreLine(1) = OrigPreLine1
!
	return
	end subroutine GiveWarning

	subroutine PrePostErrorWarning(dim, CharArray)
! Input
	integer :: dim
	character (len=*) :: CharArray(dim)
! Local
	integer i

	do i=1,dim
	  if(CharArray(i).eq.' ') return
	  write(*,fmt) trim(CharArray(i))
	end do ! i

	return
	end subroutine PrePostErrorWarning
 

	subroutine GetElasticityInfo( ThisFunction, SElastPresent, TElastPresent, &
		SElastField, TElastField, OtherFieldsPresent)
! KP added this for use with GEMPACK conversion
! Input
	type (function), pointer :: ThisFunction
! Output
! OtherFieldsPresent tells if are any fields besides s: and t: fields
	logical :: SElastPresent, TElastPresent, OtherFieldsPresent
! Each of these fields is nul if the elasticity is not present
	type (field), pointer :: SElastField, TElastField
! Local
	type (field), pointer :: FirstField, ThisField
!.NO	integer :: NumField
!.NO	logical :: Continue
	character (len=512) ElastText

	SElastPresent = .false.
	TElastPresent = .false.
	nullify(SElastField)
	nullify(TElastField)
	OtherFieldsPresent = .false.
!
	FirstField => ThisFunction%firstrecord%firstfield

	ThisField => FirstField%next
!.NO	NumField = 0
!.NO	Continue = .true.

! s: and t: must be found in the first two fields
!.NO	do while(associated(ThisField).and.(NumField.le.2).and.Continue)
	do while(associated(ThisField))
!.NO	  NumField = NumField + 1
	  if(ThisField%label.eq.'s') then
	    SElastPresent = .true.
	    SElastField => ThisField
	  else if(ThisField%label.eq.'t') then
	    TElastPresent = .true.
	    TElastField => ThisField
	  else
! Have bypassed places where s: and t: might be found
!.NO	    Continue = .false.
	    OtherFieldsPresent = .true.
	  end if
	    
	  ThisField => ThisField%next
	end do	

	end subroutine GetElasticityInfo
	
! Next added by KP to add labelinfo to symbols
	subroutine FindLabelFromLine(ThisLine, ThisLabelInfo)
! Input
	type (line), pointer :: ThisLine
! Output
	character (len=*), intent(out) :: ThisLabelInfo
! Local
	integer :: BangPos
!
	BangPos = index(ThisLine%text, '! ')
	if (BangPos > 0) then
!  Omit ! and following space
	  ThisLabelInfo = ThisLine%text(Bangpos+2:)
	else
	  ThisLabelInfo = ' '
	end if

	return
	end subroutine FindLabelFromLine

	subroutine InitializeMessages(Which, NumToSkip)  ! KP
! Input
	character (len=*), intent(in) :: Which  ! 'all' (all), 'pre' (pre) or 'post' (post)
! Only initialize those from position NumToSkip+1
	integer, intent(in) :: NumToSkip
! Local
	integer :: i

	if(Which.eq.'all'.or.Which.eq.'pre') then
	  do i=NumToSkip+1,SIZE(Messages%PreLine) 
	    Messages%PreLine(i) = ' '
	  end do
	end if
	if(Which.eq.'all'.or.Which.eq.'post') then
	  do i=NumToSKip+1,SIZE(Messages%PostLine)
	    Messages%PostLine(i) = ' '
	  end do
	end if
          
	return
	end subroutine InitializeMessages

	logical function EqualStringsIgnoreCase(str1, str2)  ! KP
! Input
	character (len=*), intent(in) :: str1, str2
! Local copies
	character (len=255) :: TempStr1, TempStr2
!
	TempStr1 = trim(str1)
	TempStr2 = trim(str2)
	call lower(TempStr1)
	call lower(TempStr2)

	EqualStringsIgnoreCase = TempStr1.eq.TempStr2

	return
	end function EqualStringsIgnoreCase

	character (len=10) function PosintToString(Posint) ! KP
! Converts a positive integer to a string
! Input
	integer, intent(in) :: Posint
! Local
	character (len=10) :: str
!
	if(posint.le.0) then
	  call InternalErrorStop('Integer is not positive')
	else if(posint.lt.10) then
	  write(str, '(I1)') Posint
	else if(posint.lt.100) then
	  write(str, '(I2)') Posint
	else if(posint.lt.1000) then
	  write(str, '(I3)') Posint
	else if(posint.lt.10000) then
	  write(str, '(I4)') Posint
	else if(posint.lt.100000) then
	  write(str, '(I5)') Posint
	else if(posint.lt.1000000) then
	  write(str, '(I6)') Posint
	else if(posint.lt.10000000) then
	  write(str, '(I7)') Posint
	else if(posint.lt.100000000) then
	  write(str, '(I8)') Posint
	else if(posint.lt.1000000000) then
	  write(str, '(I9)') Posint
	else if(posint.lt.10000000000) then
	  write(str, '(I10)') Posint
	else
	  call InternalErrorStop('Integer is too large. Routine PosintToString')
	end if

	PosintToString = str

	return
	end function PosintToString

	subroutine ProcessCommandLine(CMDLine)
	character (len=*) :: CMDLine

	integer :: ic, jc

!	Command line syntax:
!
!			mgemain -gms <inputfile>
!		or	mgemain -gp  <inputfile>
!		or	mgemain -mge <inputfile>

	call lower(CMDLine)

	ConvertToGEMPACK = (index(CMDLine,'-gp' )  > 0)
	ConvertToGAMS    = (index(CMDLine,'-gams') > 0)
	ConvertToMGE     = (index(CMDLine,'-mge')  > 0)

	call InitializeMessages('all',0)
	Messages%PostLine(1) = 'This occurred while processing the command line.'
	if (	ConvertToGEMPACK.and.(ConvertToGAMS   .or.ConvertToMGE) .or. &
		ConvertToGAMS   .and.(ConvertToGEMPACK.or.ConvertToMGE) .or. &
		ConvertToMGE    .and.(ConvertToGAMS   .or.ConvertToGEMPACK) ) &
		call ErrorStop('Multiple targets specifed -- at most one of "-gp" or "-gams" or "-mge" ')

! [07/04/04] Set ConvertToGEMPACK true by default ! KP
	if(.not.ConvertToGAMS .and. .not.ConvertToMGE) then
	  ConvertToGEMPACK = .true.
	end if

	ic = 0
	if (ConvertToGEMPACK) ic = index(CMDLine,'-gp')
	if (ConvertToGAMS)    ic = index(CMDLine,'-gams')
	if (ConvertToMGE)     ic = index(CMDLine,'-mge')
	if (ic.gt.0) then
	  jc = ic + index(CMDLine(ic:),' ') - 1
	  CMDLine (ic:jc) = ' '
	endif

	InputFileName = adjustl(CMDLine)

	return

	end subroutine ProcessCommandLine


	integer function GetArgDim(ThisArg)  !! Counts number of sets in ThisArg  ! Added by KP
!Input
	type(domain), pointer :: ThisArg
! local
	type(domain), pointer :: CurrentArg
	type(symbol), pointer :: CurrentSet

	CurrentArg => ThisArg
	GetArgDim = 0
	do while(associated(CurrentArg))
	  GetArgDim = GetArgDim + 1

	  CurrentSet => CurrentArg%s
!	  write(*,*)'GetArgDim. CurrentSet%name is:', CurrentSet%name

	  CurrentArg => CurrentArg%next
	end do ! while(associated(CurrentArg))
	
	return
	end function GetArgDim

        subroutine indices
        type (function), pointer :: fn
        type (field), pointer :: fnf
        type (nest), pointer :: n
        type (netputlist), pointer :: e
	type (endowmentlist), pointer :: el
        integer :: i
        logical :: outputs

!       Assign indices for all nests and netputs in $prod blocks and all
!       endowments in $demand blocks:

        i = 0
        fn => functionlist
        do while (associated(fn))
          fnf => fn%firstrecord%firstfield
          if (fnf%label.ne.'$prod' .and. fnf%label.ne.'$demand') then
            fn => fn%next
            cycle
          end if
	  if (fnf%label.eq.'$prod') then
            n => fn%inputs
            outputs = .false.
            do while (associated(n))
              i = i + 1
              n%index = i
              e => n%elements
              do while (associated(e))
                i = i + 1
                e%np%index = i
                e => e%next
              end do
              call nextnest(n)
              if (.not.associated(n)) then
                if (outputs) exit
                outputs = .true.
                n => fn%outputs
              end if
            end do
	  else
            el => fn%endowments
            do while (associated(el))
              i = i + 1
              el%e%index = i
              el => el%next
            end do
          endif
          fn => fn%next
        end do
        return
        end subroutine indices

	subroutine rewritemodel
	type (symbol), pointer :: s
        type(function), pointer :: fn
        type(field), pointer   :: f, fnf
	type (nest), pointer   :: n
	type (nestlist), pointer   :: nl
	type (record), pointer :: r
	type (netputlist), pointer :: npl
	character (len=1) :: prefix
	character (len=255) :: price, descript, rec
	character (len=1024) :: quant
	logical :: outputs
	integer :: i, k, ic, jc
        character (len=9) :: vtype(4)=(/'sector   ','commodity','consumer ','auxiliary'/)
	integer :: nsymbol(4)=0

!	First, load all the declared variables:

        do k=1,4
          s => idlist
          do while (associated(s))
            if (s%type.eq.vtype(k)) then
	      nsymbol(k) = nsymbol(k) + 1
              rec = trim(adjustl(s%lines%line%text))
              ic = index(rec,'!')
              if (ic.eq.0) then
	        call writegms(vtype(k),rec,' ')
              else
	        call writegms(vtype(k),rec(1:ic-1),rec(ic:))
              end if
            end if
            s => s%next
          end do
        end do

!	Then load price and quantity variables for all the nests:

	fn => functionlist
	do while (associated(fn))
	  fnf => fn%firstrecord%firstfield
	  if (fnf%label.ne.'$prod') then
	    fn => fn%next
	    cycle
	  end if
	  outputs = .false.
	  n => fn%inputs
	  do while (associated(n))
	    write(*,fmt) 'Nest check: '//trim(n%ident)
	    write(*,*) associated(n%parent)
	    write(*,*) associated(n%subnests)
	    write(*,*) associated(n%nlist)
	    n => n%next
	  end do
	  n => fn%inputs
	  do while (associated(n))
	    if (associated(n%parent)) then
	      call writeparam('p'//'_',n%index,fnf%arg,price)
	      if (outputs) then
	        call writeparam('s'//'_',n%index,fnf%arg,quant)
	        descript = ' ! Supply index for output nest '//trim(n%ident)//' from sector '//trim(fnf%data)
		call writegms('sector',quant,descript)
		descript = ' ! Price index for output nest '//trim(n%ident)//' from sector '//trim(fnf%data)
		call writegms('commodity',price,descript)
	      else
	        call writeparam('d'//'_',n%index,fnf%arg,quant)
	        descript = ' ! Demand index for input nest '//trim(n%ident)//' in sector '//trim(fnf%data)
		call writegms('sector',quant,descript)
		descript = ' ! Price index for input nest '//trim(n%ident)//' in sector '//trim(fnf%data)
		call writegms('commodity',price,descript)
	      end if
	    end if
	    write(*,*) ' Nest: '//trim(n%ident)
	    call nextnest(n)
	    if (.not.associated(n) .and. .not.outputs) then
	      n => fn%outputs
	      outputs = .true.
	    end if
	  end do
	  fn => fn%next
        end do

	write(iomge,fmt) '$model:'//trim(modelname)
	if (nsymbol(1).gt.0) then
	  write(iomge,'(/a)') '$sectors:'
	  call writetype(iomge,8,'sector')
	end if
	write(iomge,'(/a)') '$commodities:'
	call writetype(iomge,8,'commodity')
	write(iomge,'(/a)') '$consumers:'
	call writetype(iomge,8,'consumer')
	if (nsymbol(4).gt.0) then
	  write(iomge,'(/a)') '$auxiliary:'
	  call writetype(iomge,8,'auxiliary')
	end if

!	Then write out all the demand functions:

	fn => functionlist
	do while (associated(fn))
	  r => fn%firstrecord
	  fnf => r%firstfield
	  if (fnf%label.ne.'$demand') then
	    fn => fn%next
	    cycle
	  end if
	  write(iomge,*)
	  do i=1,fn%nrecord
	    write(iomge,fmt) trim(r%line%text)
	    r => r%next
	  end do
	  fn => fn%next
	end do

!	Then write out all the production functions:

	fn => functionlist
	do while (associated(fn))
	  r => fn%firstrecord
	  fnf => r%firstfield
	  if (fnf%label.ne.'$prod') then
	    fn => fn%next
	    cycle
	  end if
	  write(iomge,*)
	  write(iomge,fmt) trim(fnf%label)//':'//trim(fnf%data)//' s:'//trim(fn%inputs%value)//' t:'//trim(fn%outputs%value)

	  npl => fn%outputs%elements
	  do while (associated(npl))
	    write(iomge,'(5x,a)') trim(adjustl(npl%np%f%line%text))
	    npl => npl%next
	  end do
	  nl => fn%outputs%subnests
	  do while (associated(nl))
            call writeparam('p_',nl%n%index,fnf%arg,price)
            call writeparam('r0_',nl%n%index,fnf%arg,quant)
	    write(iomge,'(5x,a)') 'i:'//trim(price)//' q:'//trim(quant)
	    nl => nl%next
	  end do

	  npl => fn%inputs%elements
	  do while (associated(npl))
	    write(iomge,'(5x,a)') trim(adjustl(npl%np%f%line%text))
	    npl => npl%next
	  end do
	  nl => fn%inputs%subnests
	  do while (associated(nl))
            call writeparam('p_',nl%n%index,fnf%arg,price)
            call writeparam('c0_',nl%n%index,fnf%arg,quant)
	    call substitute(quant)
	    write(iomge,'(5x,a)') 'i:'//trim(price)//' q:'//trim(quant)
	    nl => nl%next
	  end do

	  n => fn%inputs
	  outputs = .false.
	  do while (associated(n))
	    if (associated(n%parent)) then
	      call writeparam('p'//'_',n%index,fnf%arg,price)
	      if (outputs) then
	        call writeparam('s_',n%index,fnf%arg,quant)
		write(iomge,'(/a)') '$prod:'//trim(quant)//' t:'//trim(n%value)
	        call writeparam('r0_',n%index,fnf%arg,quant)
		call substitute(quant)
		write(iomge,'(5x,a)') 'i:'//trim(price)//'  q:'//trim(quant)
	      else
	        call writeparam('d_',n%index,fnf%arg,quant)
		write(iomge,'(/a)') '$prod:'//trim(quant)//' s:'//trim(n%value)
	        call writeparam('c0_',n%index,fnf%arg,quant)
		call substitute(quant)
		write(iomge,'(5x,a)') 'o:'//trim(price)//'  q:'//trim(quant)
	      end if
	      npl => n%elements
	      do while (associated(npl))
	        rec = trim(adjustl(npl%np%f%line%text))
		ic = index(rec,trim(n%ident)//':')
		if (ic.eq.0) stop 'Cannot find nest identifier text?'
		jc = ic+len_trim(trim(n%ident)//':')
		write(iomge,'(5x,a)') rec(1:ic-1)//rec(jc:)
		npl => npl%next
	      end do
	      nl => n%subnests
	      do while (associated(nl))
	        call writeparam('p_',nl%n%index,fnf%arg,price)
		call writeparam('c0_',nl%n%index,fnf%arg,quant)
		call substitute(quant)
		write(iomge,'(5x,a)') 'i:'//trim(price)//' q:'//trim(quant)
		nl => nl%next
	      end do
	    end if
	    call nextnest(n)
	    if (.not.associated(n) .and. .not.outputs) then
	      n => fn%outputs
	      outputs = .true.
	    end if
	  end do

	  fn => fn%next
	end do

	return
	end subroutine rewritemodel

	subroutine substitute(item)
	character (len=*) :: item
	type (gmsrecord), pointer :: r	
	integer ic, jc, kc, lc

	if (item(1:1).ne.'(') item = '('//trim(item)//')'
	write(*,fmt) 'Substitute:'//trim(item)
	r => gmsrecords
	do while (associated(r))
	  if (r%type.eq.'costdef') then
	    write(*,fmt) 'Checking: '//trim(r%left)
	    ic = index(r%left,'=')
	    if (ic.eq.0) stop 'Error : no sign found.'
	    jc = index(item,trim(r%left(1:ic-1)))
	    if (jc.gt.0) then
	      write(*,fmt) 'Substituting'
	      kc = jc + len_trim(r%left(1:ic-1))
	      lc = index(r%left,';') - 1
	      if (jc.gt.1) then
	        item = item(1:jc-1)//trim(adjustl(r%left(ic+1:lc)))//item(kc:)
	      else
	        item = trim(adjustl(r%left(ic+1:lc)))//item(kc:)
	      end if
	      r => gmsrecords
	    else
	      r => r%next
	    end if
	  else
	    r => r%next
	  end if
	end do
	return
	end subroutine substitute

        subroutine nextnest(n)
        type (nest), pointer :: n

        if (associated(n%subnests)) then
          n => n%subnests%n
          return
        endif

!       Keep looking until we have found another the 
!       next nest:

        do 
          if (associated(n%nlist)) then
            if (associated(n%nlist%next)) then
              n => n%nlist%next%n
              return
            else

!       No more items in the list to be processed, so we
!       shift focus to the parent nest and continue looking:

              n => n%parent
            end if
          else

!       No parent list found which means that we have completed
!       our traversal of the entire tree:

            nullify(n)
            return
          end if
        end do

        end subroutine nextnest

	subroutine writegms(rectype,lrec,rrec)
	character (len=*) :: rectype, lrec, rrec
	type (gmsrecord), pointer :: r
	allocate(r)
	r%type = rectype
	r%left = adjustl(lrec)
	r%right = adjustl(rrec)
	nullify(r%next)
	if (.not.associated(gmsrecords)) then
	  gmsrecords => r
	  lastrecord => r
	else
	  lastrecord%next => r
	  lastrecord => r
	end if
	nullify(r)
	return
	end subroutine writegms

	subroutine writetype(ioout, tabin,type)
	integer :: ioout
	character (len=*) :: type
	type (gmsrecord), pointer :: r	
	character (len=32) :: fmt
	integer :: nc, tabin
	nc = 0
	r => gmsrecords
	do while (associated(r))
	  if (r%type.eq.type)  nc = max(nc, len_trim(r%left))
	  r => r%next
	end do
	nc = nc + tabin
	if (tabin.eq.0) then
	  write(fmt,'(a,i0,a)') '(a,t',nc+2,',a)'
	else
	  write(fmt,'(a,i0,a,i0,a)') '(',tabin,'x,a,t',nc+2,',a)'
	endif
	r => gmsrecords
	do while (associated(r))
	  if (r%type.eq.type) write(ioout,fmt) trim(r%left), trim(r%right)
	  r => r%next
	end do
	return
	end subroutine writetype
	

!----------------------------------------------------------------------
! [06/04/04] Old subroutine ProcessCommandLine as written by KP.
! Not used now. But may want some of these features later. So
!  change its name and note that it is not called at present.
	subroutine ProcessCommandLine_OldVersion(CMDLine)
! User may specify COnvertToGEMPACK or ConvertToGAMS on command line
! User may specify name InputFileName of input file on command line
! Input
	character (len=*), intent(in) :: CMDLine
! Local
	integer :: StartPos, TrueStart, TrueEnd, llen
	integer :: NumWordsFound
	logical :: Found, Finished, Answer1, Answer2
	
! Initialize things
	ConvertToGEMPACK = .false.
	ConvertToGAMS = .false.
	InputFileName = ' '

	llen=LEN_TRIM(CMDLine)

! Go through CMDLine checking what if anything has been set.
!  Need to pull off words one at a time. Quotes "" can enclose spaces in a word
	StartPos=1
	Finished = .false.
	NumWordsFound = 0

	do while(StartPos.le.llen.and. (.not.Finished))
	  call GetNextWord(CMDLine, StartPos, &
	    found, TrueStart, TrueEnd)
	  if(.not.found) then
	    Finished = .true.
	  else
! Here have found another word
! Process it.
	    NumWordsFound = NumWordsFound + 1

!	    write(*,*) 'Word number ', NumWordsFound, 'is: ', CMDLine(TrueStart:TrueEnd)

! [31/03/04] For the present, always ConvertToGEMPACK. So ignore this "word 1" case"
	    ConvertToGEMPACK = .true.

!LATER	    if(NumWordsFound.eq.1) then
! Expect "-gp" or "-gams"
!LATER	      Answer1 = EqualStringsIgnoreCase(CMDLine(TrueStart:TrueEnd), '-gp')
!LATER	      Answer2 = EqualStringsIgnoreCase(CMDLine(TrueStart:TrueEnd), '-gams')
!LATER	      if(Answer1) then
!LATER		ConvertToGEMPACK = .true.
!LATER	      else if(Answer2) then
!LATER		ConvertToGAMS = .true.
!LATER	      else
! Unexpected word
!LATER		call InitializeMessages('all',0)
!LATER		Messages%PostLine(1) = 'This occurred while processing the command line.'
!LATER		call ErrorStop('Unexpected first word on command line. Expect "-gp" or "-gams".')
!LATER	      end if
!LATER	    else if(NumWordsFound.eq.2) then

! [31/03/04] Just process one word at present
	    if(NumWordsFound.eq.1) then
! Expect input file name
	      InputFileName = CMDLine(TrueStart:TrueEnd)
	    end if
	  end if

!LATER At present, only process at most 2 words
!LATER	  if(NumWordsFound.ge.2) then
! [31/03/04] Just process one word at present
	  if(NumWordsFound.ge.1) then
	    Finished = .true.
	  end if

	  if(Found) then
	    StartPos=TrueEnd+1
	  end if

	end do ! while(StartPos.le.llen.and. (.not.Finished))

	return
	end subroutine ProcessCommandLine_OldVersion

! [06/04/04] Subroutine GetNextWord was called by old version of
!  routine ProcessCommandLine as written by KP.
! This subroutine GetNextWord is not called at present.
!  [But may want some of these features later.]
	subroutine GetNextWord( &
! GetNextWord from a string (eg the CommandLine)
! This is based on GEMPACK routine GPGNW3.
! Quotes " can be used to start and end a word which may contain spaces 
! Input
		Line, StartPos, &
! Output
	  	Found, TrueStart, TrueEnd)
! Input	
	character (len=*), intent(in) :: Line
	integer, intent(in) :: StartPos
! Output
	logical, intent(out) :: Found
! If found, these indicate start and end of next word
	integer, intent(out) :: TrueStart, TrueEnd 	


! Local
	character (len=512) line2
	character :: ch
	integer :: i, qupos, n1, llen
	logical :: quotok
	logical :: Finished, FoundStart

	Finished = .false.
	FoundStart = .false.

	llen=LEN_TRIM(LIne)

! Quick test in case rest of line is all blanks
!
	if(StartPos.le.llen) then
	  if(line(StartPos:).eq.' ') then
	    found=.FALSE.
	    Finished = .true.
!GP	      Return
!GP	    goto 400
	  end if
	end if

	if(.not.Finished) then
	  do I = StartPos, llen
	    ch=line(i:i)
! Quotes
	    if(.not.FoundStart.and.ch.eq.'"') then
!		Word starts with " and ends with next " on this line
!		Then TrueStart, TrueEnd point to the quotes.
!
	      qupos=i
	      quotok=.TRUE.
!
!		QUOTOK is set false if is not matching " later on this line
!
	      if(qupos.eq.llen) then
		quotok=.FALSE.
	      else
	        TrueStart=qupos
	        n1=INDEX(line(qupos+1:),'"')
	        if(n1.gt.0) then
		  TrueEnd=qupos+n1
	        else
		  quotok=.FALSE.
	        end if
	      end if
!
	      if(quotok) then
!
		found=.TRUE.
!GP		goto 400
		FoundStart = .true.
!
	      else
!
!		Here no matching ending quote.
!		Message
!
		line2=' '
		line2(qupos+1:)='?'

!GP		call GPTO2('\n%t','\n%t',line,line2,0,0,errmes)
!GP	        if(errmes(1:1).eq.'E') goto 15000
!GP		call GPTO1(
!GP     * '\n \%\% No matching quote for quote in position %i.',
!GP     *		  ' ', qupos, errmes)
!GP	        if(errmes(1:1).eq.'E') goto 15000

		call InitializeMessages('all',0)
		Messages%PostLine(1) = 'This occurred while processing the command line.'
		call ErrorStop('Unmatched quotes "" on command line')
!
		found=.false.
!GP		goto 400
		Finished=.true.
!
	      end if
!
	    end if ! if(.not.FoundStart.and.ch.eq.'"') then

! TABs ??   if(ch.ne.' '.and.ch.ne.tab) then
	    if(.not.Finished.and..not.FoundStart.and.ch.ne.' ') then
	      TrueStart=i
	      FoundStart = .true.
! GP	      goto 200
	    end if
	  end do  ! do i

	end if ! if(.not.Finished) then

	if(.not.Finished.and.FoundStart) then
!
!       **** Here we have found the start of a new substring. ****
!
          FOUND=.TRUE.
!
!       **** Find the end of this string. Set TEND equal to this. ****
!
          DO I = TrueStart,llen
	    ch=line(i:i)
!
! TABs??  if(ch.eq.' '.or.ch.eq.tab) then
	    if(ch.eq.' ') then
	      TrueEnd=i-1
	      Finished=.true.
	      Found = .true.
!GP	      goto 400
	    end if
	  end do ! do i
!
!       **** if this point reached end of string is at end of input line. ****
!
	  if(.not.Finished) then
	    TrueEnd = llen
	      Finished=.true.
	      Found = .true.
	  end if

	end if ! if(.not.Finished.and.FoundStart) then

	return
	end subroutine GetNextWord
! This subroutine GetNextWord is not called at present.
!----------------------------------------------------------------------

end module mpsge

